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Abstract 


Acoustic  landmine  detection  is  accomplished  by  using  a  loud  speaker  to  generate  airborne  source 
low-frequency  waves  that  are  transmitted  to  the  soil  above  a  buried  landmine  target.  At  a  specific 
frequency,  the  landmine  will  “vibrate”  at  resonance,  imparting  an  enhanced  velocity  on  the  soil 
particles  above  it  at  the  surface  that  is  detected  by  a  scanning  Laser  Doppler  Vibrometer  system. 
If  the  soil  surface  velocity  profiles  measured  from  these  experiments  could  be  predicted  mathe¬ 
matically  under  a  variety  of  conditions,  the  physical  system  would  be  able  to  accurately  detect 
landmines  in  more  challenging  environments. 

The  mathematical  modeling  of  the  buried  landmine  detection  problem  involved  wave  propaga¬ 
tion  in  a  layered  waveguide  in  the  presence  and  absence  of  a  buried  circular  target.  In  this  study, 
emphasis  was  placed  on  acoustic  to  seismic  coupling  of  an  airborne  continuous  wave  point  source 
into  the  soil.  Soil  resonances  were  calculated  with  a  model  that  represents  the  soil  as  a  finite,  fluid- 
hlled  rigid  porous  layer  below  a  finite  atmospheric  layer.  This  two-layer  waveguide  incorporated 
density  and  sound  speed  in  both  the  soil  and  atmosphere,  which  was  adjusted  based  on  soil  type, 
compactness,  and  moisture  content  in  both  the  air  and  soil.  An  analytic  solution  of  the  two-layer 
waveguide  problem  involved  solving  the  Helmholtz  equation  in  cylindrical  coordinates  in  both  lay¬ 
ers  along  with  using  a  delta  function  point  source  to  simulate  a  compact  loudspeaker  in  the  upper 
layer.  Boundary  conditions  along  with  conditions  of  orthogonality  were  used  to  obtain  complicated 
analytical  expressions  for  the  eigenvalues  and  eigenfunctions  of  the  problem.  A  MATLAB™  pro¬ 
gram  was  used  to  numerically  solve  for  the  eigenvalues  and  plot  solutions  of  pressure  and  particle 
velocity  vs.  frequency  and  spatial  variables. 

In  the  presence  of  the  buried  circular  target,  membrane  resonances  were  predicted  by  similar 
modeling  techniques.  The  top  plate  of  the  buried  landmine  was  modeled  as  a  circular  elastic 
membrane  stretched  flush  over  a  cylindrical  cavity  in  a  rigid  substrate  beneath  the  porous  layer. 
The  Helmholtz  equation  was  again  solved  in  the  atmospheric  layer  using  cylindrical  coordinates 
and  a  point  source.  The  homogeneous  Helmholtz  equation  was  again  used  in  the  porous  layer,  with 
a  Green’s  function  representation  of  the  membrane  response. 

In  the  soil  resonance  predictions,  pressure  was  plotted  as  a  function  of  frequency,  and  the 
resonances  appear  as  local  minimums  and  maximums.  A  MATLAB™  user  interface  was  created 
to  allow  researchers  to  access  the  soil  resonance  information  particular  to  their  experiment.  The 
resonances  of  the  membrane-soil  model  were  determined  using  MATLAB™;  describing  the  effects 
of  frequency,  depth,  density,  and  sound,  to  include  radius  and  elastic  parameters  of  the  membrane. 
Comparison  of  the  results  (involving  the  fluid  surface  particle  velocity  profiles  across  the  target) 
has  been  made  with  experiments  reported  in  the  literature  to  evaluate  the  model’s  usefulness. 
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Chapter  2 
Introduction 


Landmines  first  appeared  in  warfare  in  the  third  century  BC,  however,  the  first  use  of  non-metallic 
landmines  was  in  World  War  I.  Due  to  metal  shortages  at  the  end  of  the  war,  Germany  began 
implementing  primitive  wooden  antitank  mines  against  Allied  forces.  However,  the  Gennans 
reverted  to  metallic  mines  by  the  start  of  the  Second  World  War.  By  1942,  Allied  landmine  de¬ 
tection  technology  had  advanced  significantly:  from  a  long  stick  with  a  metal  point  used  to  prod 
the  ground,  to  the  conventional  metal  detectors  common  today.  The  natural  solution  was  to  begin 
making  landmines  from  other  materials,  notably  wood,  Bakelite  (a  simple  plastic),  and  glass.  For 
example,  the  Gennan  Topfmine  was  an  antitank  mine  made  of  plastic,  wood,  glass,  and  cardboard; 
it  required  over  300  pounds  of  pressure  to  be  applied  to  its  top  to  discharge  the  mine.  Approx¬ 
imately  800,000  of  these  mines  were  constructed  during  the  last  year  of  the  war.1  Similar  to  the 
Topfmine,  the  Glasmine-43  was  the  non-metallic  antipersonnel  landmine  developed  by  the  Ger¬ 
mans.  Each  mine  was  laid  precisely  and  strategically,  and  the  Germans  kept  meticulous  records  of 
their  minefields  throughout  the  war.  Non-metallic  landmines  were  also  used  in  both  the  Korean 
and  Vietnam  Wars.  Frequently,  the  landmines  that  were  laid  in  both  of  these  countries  and  around 
the  world  during  the  Cold  War  era  were  not  documented,  and  the  rare  instances  of  documentation 
are  no  longer  accurate  due  to  soil  shifting  over  time. 

Today,  Afghanistan  is  one  of  the  most  heavily  mined  nations  in  the  world.2  The  first  wide¬ 
spread,  massive  deployment  of  landmines  in  Afganistan  was  on  the  Pakistan  and  Iran  borders  by 
the  Soviet  Anny  to  enforce  a  policy  of  area-denial  around  1979.  U.S.  soldiers  continue  to  lose  their 
lives  to  landmines  in  Operation  Enduring  Freedom.3  Landmines  also  affect  the  civilian  population. 
Red  Cross  estimates  indicate  17  in  1,000  children  in  Afghanistan  have  been  injured  or  killed  by  a 
landmine.4  The  area  that  could  not  be  controlled  by  the  Soviets  in  1979  is  still  unavailable  to  the 
people  of  Afghanistan,  who  are  mostly  fanners  in  need  of  more  crop  land.  A  cost  effective  method 
for  landmine  detection  and  removal  is  essential  for  political  and  economic  reasons  in  this  nation. 

Landmines  are  used  in  a  wide  variety  of  political  and  economic  situations,  including  wars,  bor- 

'[28]  pg.  113. 
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der  disputes,  conflicts,  and  coup  d’etats.5  The  broad  usage  of  landmines  is  partially  made  possible 
by  a  variety  of  methods  of  detonation,  including  a  change  in  pressure  from  personnel  or  tanks 
transiting  over  them,  trip  wires,  or  remote  firing  from  a  predetermined  location.6  Additionally,  the 
low  cost  (between  $3-5  U.S.)  7  and  effectiveness  of  killing  or  maiming  humans  and  destroying 
machinery  has  made  the  landmine  a  tactical  military  asset.  For  factions  seeking  to  overthrow  a 
government,  the  laying  of  landmines  instills  both  chaos  and  panic  within  the  civilian  population 
that  the  government  is  seemingly  powerless  to  stop.  Landmines  are  present  in  approximately 
one-third  of  developing  nations  and  continue  to  kill  and  affect  lives  long  after  wars  and  conflicts 
have  ended.  According  to  the  United  Nations  in  1995,  over  100  million  landmines  were  in  use  in 
conflicts  throughout  the  world.8  Estimates  of  the  cost  to  remove  landmines  range  from  $200  (U.S.) 
to  over  $1000  (U.S.)  per  mine,9  making  removal  by  many  third  world  governments  economically 
impossible. 

The  acoustic  landmine  detection  problem  has  gained  much  interest  in  the  scientific  research 
community  in  landmine  detection  and  removal  over  the  last  ten  years.  Some  landmine  detection 
equipment  has  not  been  updated  significantly  since  World  War  II,  and  the  current  techniques  us¬ 
ing  ground  penetrating  radar  still  need  to  be  perfected  for  finding  non-metallic  landmines  due  to 
the  high  rate  of  false  alarms.  Radar  also  generates  large  numbers  of  false  alarms  due  to  soil 
inhomogenities,  that  is,  it  cannot  discriminate  between  the  frequency  response  of  the  soil  inho- 
mogenities  and  landmine  resonances.  To  further  complicate  matters,  the  technology  for  detecting 
non-metallic  landmines  is  very  limited;  the  only  effective  known  methods  are  seismic  acoustic 
wave  penetration  and  vapor  detection.  Vapor  detection,  a  chemical  method  for  detection  which  an¬ 
alyzes  gas  emissions  from  the  mine,  is  slow,  expensive,  and  dangerous.  Additionally,  this  method 
only  indicates  that  a  mine  was  present  at  some  time  in  the  soil,  and  cannot  predict  its  exact  loca¬ 
tion,  making  it  particularly  dangerous.  Tests  require  new  chemicals  that  are  not  reusable,  which 
increase  the  cost  of  landmine  detection. 

Sabatier  and  Waxler’s  work  at  the  National  Center  for  Physical  Acoustics  (NCPA),  supported 
by  the  U.S.  Army,  Navy,  and  Marine  Corps  de-mining  programs,  suggests  that  acoustic  and  seis¬ 
mic  technology  may  be  the  most  cost  effective  and  accurate  way  to  detect  landmines.10  NCPA’s 
work  on  acoustic  to  seismic  landmine  detection  is  one  of  the  most  promising  methods  for  detecting 
non-metallic  landmines  at  shallow  burial  depths  with  low  false-alarm  rates.  A  loud  speaker  acts 
as  an  airborne  source,  generating  low-frequency  waves  which  enter  the  soil  over  a  buried  target. 
At  a  specific  frequency  known  as  the  resonant  frequency,  the  mine  will  “vibrate”  with  a  relatively 
large  amplitude,  imparting  a  certain  velocity  on  the  soil  particles  above  it  that  is  detected  at  the 
surface  by  a  scanning  Laser  Doppler  Vibrometer  system.  In  field  tests  with  the  Valmara  VS  1.6 
plastic-cased,  anti-tank  landmine,  the  system  has  worked  very  well,  with  a  probability  of  detection 
exceeding  95%  in  desert  environments.11  There  was  only  one  false  alarm  during  the  experiment, 


5[7],pg. 

126. 

6[21],  Pg- 

18. 

7[14],pg. 

1. 

8[14],  pg. 

1. 

9[14],  pg. 

1. 

10[19],  pg. 

149. 

"[19],  pg. 

28. 

CHAPTER  2.  INTRODUCTION 


10 


despite  rocky  soil  in  an  Arizona  desert  test  site.  If  the  soil  surface  velocity  plots  could  be  predicted 
mathematically  under  a  variety  of  conditions,  the  physical  system  would  be  able  to  better  discrim¬ 
inate  between  the  frequency  response  of  the  false  alarms  and  the  landmine  in  more  challenging 
environments. 

Sabatier  suggests  that  the  success  of  airborne  acoustics  in  landmine  detection  depends  on  three 
assumptions.  The  first  assumption  is  that  landmines  are  significantly  more  acoustically  compli¬ 
ant  than  soil  and  other  natural  materials  (rocks,  roots,  etc.).  Second,  the  acoustic  properties  of  a 
landmine  case  are  assumed  to  further  differentiate  the  mine  through  contrast  with  the  porous  soil. 
Last,  the  interface  between  the  soil  and  mine  is  thought  to  be  non-linear,  which  implies  that  the 
interface  relation  is  complicated,  further  contrasting  the  mine  from  the  soil  (although  it  is  assumed 
linear  in  this  paper  for  simplicity).  The  nonlinearity  of  soil  vibration  and  soil  interacting  with  the 
top  plate  of  the  mine  case  is  a  future  problem  that  could  be  developed  from  this  Trident  Research 
Project  down  the  road.  The  use  of  algorithms  implemented  in  the  experimental  testing  has  resulted 
in  "almost  perfect  detection  and  zero  false  alann  rates."12 

Sabatier  also  indicates  that  acoustic  to  seismic  mine  detection  works  best  in  a  dry,  sandy  soil 
environment.13  However,  he  also  states  that  landmine  detection  is  immune  to  atmospheric  moisture, 
weather,  other  acoustic  (noise)  sources,  and  natural  materials,  according  to  his  experiments.14  Any 
man-made  objects  (soda  cans,  etc.)  scattered  in  the  soil  may  be  detected  as  false  positives.  Large 
amounts  of  vegetation  or  grass,  as  well  as  frozen  soil,  appear  to  compromise  detection,  since  the 
majority  of  acoustic  energy  is  absorbed  or  deflected.  Despite  these  limitations,  Sabatier  is  highly 
encouraged  by  the  results  that  he  has  been  able  to  obtain  in  such  a  short  time,  and  reports  that  the 
system  should  be  ready  for  implementation  in  4-5  years.15 

The  goal  of  this  project  is  to  theoretically  predict  the  acoustic  to  seismic  soil  vibration  velocity 
profiles  that  are  related  to  the  experiments  reported  by  Sabatier  and  others.  The  landmine  scat¬ 
tering  problem  has  not  been  solved  analytically  because  of  its  complexity,  but  has  been  solved 
computationally  for  one  special  case.16  The  two  mathematical  scattering  geometries  include  (1) 
an  atmospheric  layer  over  a  single  layer  of  soil  bounded  by  a  rigid  substrate,  and  (2)  the  simpli¬ 
fied  geometry  demonstrating  propagation  through  a  homogeneous  soil  onto  a  rigid  surface  with 
a  circular  clamped  vibrating  membrane  over  a  cavity  in  the  rigid  surface.  Secondly,  the  models 
include  soil  variations,  angle  of  incidence,  effects  of  acoustic  wave  scattering,  size  and  other  char¬ 
acteristics  of  the  target,  and  depth  of  the  target  from  the  surface  through  the  programs  devised  for 
each  geometry  (with  the  exception  of  soil  variation,  which  is  incorporated  individually  into  the 
program). 

The  mathematical  modeling  of  the  wave  equation  includes  implementation  of  various  compu¬ 
tational  tools,  including  MATLAB™  and  COMSOL™.  COMSOL™  is  a  finite  element  partial 
differential  equation  solver,  which  includes  a  specialized  acoustics  interface.  COMSOL™  was 
used  as  a  method  of  comparison  to  evaluate  the  effectiveness  of  the  model.  Both  methods  are 

12[19],  pg-  151. 

13[27],  pg-  1-2. 

14[27],  pg.  1-2. 

15[19],  pg.  151. 

16[27], 
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compared  to  verify  accuracy.  Finally,  a  user  friendly  interface  program  allows  researchers  in  the 
acoustic  field  to  use  the  mathematical  profile  predictions  in  this  project.  The  user  interface  will 
hopefully  contribute  to  future  research  in  the  acoustics  detection  of  the  buried  landmine  problem. 


Chapter  3 

Background  and  Prerequisites 


3.1  Derivation  of  the  Wave  Equation 

3.1.1  Overview 

The  derivation  of  the  wave  equation  presented  in  this  paper  follows  the  argument  in  Morse  and 
Ingard,1  and  Zwikker  and  Kosten2.  Suppose  a  fluid  infiltrates  a  rigid  porous  solid  that  resists  the 
fluid’s  flow  and  allows  for  changes  the  properties  of  the  fluid.  The  pores  of  the  solid  material  are 
randomly  interconnected,  and  it  is  assumed  that  fluid  flow  in  each  direction  has  identical  acoustic 
impedance.  Let  Tt  denote  the  porosity,  the  ratio  of  the  fluid  volume  to  the  total  volume.  Let  w 
be  the  mean  flow  velocity  for  the  porous  solid.  The  fluid  has  uniform  density  (p),  pressure  (P), 
entropy  (S),  and  temperature  (T)  in  the  absence  of  sound.  When  an  incident  sound  wave  impacts 
the  fluid,  the  pressure  changes  to  P  =  P  +  p  (x,  t),  assuming  the  wave  is  one  dimensional  and 
propagating  in  the  x  direction.  The  density,  entropy,  and  temperature  also  change  to  p  =  p+S  (x,  t), 
S  =  S  +  s(x,  t),  and  T  =  T  +  r  (x,  t ),  respectively.  Assume  that  p,  5,  s,  and  r  are  small  when 
compared  with  P,  p,  S,  and  T,  respectively.  The  only  energy  resulting  from  the  acoustic  motion 
is  mechanical,  since  other  forms  of  energy  are  eliminated  by  assuming  the  medium  is  nonviscous 
and  non-conducting. 

3.1.2  Adiabatic  Compressibility 

Due  to  the  acoustic  motion  being  completely  mechanical,  the  forces  in  the  problem  are  limited  to 
compressive  elasticity.  Morse  and  Ingard  conclude  that  at  frequencies  lower  than  1  GHz,  compres¬ 
sion  is  adiabatic.3  The  compression  is  referred  to  as  adiabatic  because  no  heat  is  transferred  to  or 
from  the  fluid  during  the  compression.  Thus,  for  an  ideal  gas,  p  =  p  (P,  S') ,  and  therefore, 

P  (P,  S)  =  p  (P,  S)  +  Ap  (p,  S)  p  +  Ap  (P,  S)  a 

‘[22],  pg.  253. 

2[29],  pg.  18. 

3 [22],  pg.  230. 
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to  the  first  order.  Since  adiabatic  compressibility  implies  s  =  0, 

P  (p,s)  =p(P,S)  +  -^pp  (P,  S)  p  =  p  +  pnsp  (3.1) 

where  the  adiabatic  compressibility  is 

1  dp 

ft S  = 

p  dP 

Since  P  —  f  ( S )  p1  for  an  ideal  gas  (where  /  (S)  is  a  function  of  entropy  and  7  is  the  ratio  of  the 

c 

specific  heat  at  constant  pressure  of  the  gas  to  the  specific  heat  at  constant  volume,  ie,  7  =  gE), 


and  therefore, 


From  equation  (3.1), 


dp 

dP 


1  =  f(S)  7 P 
1 

f(S)  7 p7”1 


7-i 


dP 

1 

^7P7'1 


_P_ 

7  P’ 


K 


s 


1 

w 


5  =  pnsp. 


(3.2) 


3.1.3  Conservation  of  Momentum 


Consider  a  small  volume  AC  =  AxAyAz  moving  in  a  vector  velocity  field  w.  The  pressure 
exerted  on  the  x-faces  of  the  cube  is 


p  (x,  y,  z,  t )  AyAz  —  p  (x  +  Arc,  y,  z,  t)  AyAz  = 


P  {x,  y,  z,t)  —  p(x  +  Ax,  y,  z,  t ) 


Ax 


AV 


d 


—p  (x,  y,  z,  t)  AV. 


The  total  pressure  exerted  on  the  entire  cube  is  then  —'VpAV.  Newton’s  second  law  gives 

dw 

pAC—  =  -VpAC, 


(3.3) 


where  V  is  the  vector  differential  operator.  Gravitational  force  is  not  included  in  this  derivation. 
(See  the  discussion  in  Appendix)  The  volume  element  follows  a  trajectory  (■ x(t),y(t),z(t )).  The 
velocity  vector  is  tangent  to  the  curve:  w  =  (wx,wy,wz)  =  (x'(t),y'(t),z'(t)).  From  the  chain 
rule, 


d_ 

dt 


w  (x (t),y(t),z(t),t) 


(iw)^  +  (|w)^  +  (iw)*'w  +  sw 

d  d  d  d 

Wx—w  +  Wy—w  +  wz— w  +  — w 

dx  dy  dz  dt 

(w  •  V)  w  +  — w 
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From  equation  (3.3) 


d 

p  (  (w  •  V)  w  +  — w  )  =  —  S/p. 
at 


Assuming  the  components  of  w  are  small,  the  linearized  equation  for  continuity  of  momentum  is 

(3.4) 


<9w 

p~m  =  -Vp 


3.1.4  Equation  of  Continuity  (or  Conservation  of  Mass) 

The  mass  in  a  small  volume  element  is 

m  =  JJJpdV, 

R 

where  m  is  mass  and  R  is  a  small  volume  element.  Therefore, 

f  -  sn%dv 


R 


dt 


%  represents  the  mass  flowing  into  or  out  of  R.  Assuming  conservation  of  mass, 

(I^r  =  -ffp*r-~ndS 

ax  dR 

where  ~n  is  the  outward  normal  vector  to  dR.  According  to  the  divergence  theorem, 

^  =  -///V  •  (pw)  dV. 


Therefore, 


and 


R 


!JJpv  =  -fjfv-(pw)dv 

Ral  R 

N  (|  +  v '  <pw>)  dv  =  °- 


Letting  R  shrink  down  to  a  point,  (x,  y,  z ),  yields 

Yi  +  V  '  (Pw)  =  °- 

Assuming  w  is  small  (and  since  5  was  already  assumed  small),  (p  +  5)  w  ~  pw.  Thus,  the  lin¬ 
earized  equation  for  conservation  of  mass  is 

dS  „  f  , 

—  +  V  •  (pw)  =  0 


and  from  equation  (3.2) 


Ks%  =  _V  ’  ('PW') ' 


(3.5) 
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3.1.5  The  One  Dimensional  Wave  Equation 

From  conservation  of  momentum,  (eqn.  (3.4)). 

dwx  dp 

^  dt  dx 

where  w  is  velocity.  The  equation  for  continuity  for  the  density  is 

d8_  _  d  ( pwx ) 
dt  dx 


(3.6) 


If  we  assume  that  the  change  in  pressure  is  related  to  the  change  in  density  by  8  (x,  t)  ~  pnsp  (x,  t) 
(eqn.  3.2),  then 

dp  dwx  dp 


PKs~r^  =  ~P  (1  +  Ksp)  -7—  -  pnswx 


dx 


(3.7) 


Since  p,  8,  and  r  are  small  in  comparison  to  P,  p,  and  T,  the  second  and  third  order  tenns  contain¬ 
ing  p,  8,  and  r  in  both  equation  (3.6)  and  equation  (3.7)  are  neglected.  Therefore,  eqns.  (3.6)  and 
(3.7)  are  now 

dp  dwx 
dx  ^  dt 


dp 

dt 


dwx 

dx 


According  to  Morse  and  Ingard,4 


The  first  of  these  states  that  a  pressure  gradient  [along  the  x-direction]  produces  an  ac¬ 
celeration  of  the  fluid;  the  second  states  that  a  velocity  gradient  [along  the  x-direction] 
produces  a  compression  of  the  fluid. 


By  differentiating  the  first  equation  by  x  and  the  second  equation  by  t  eliminates  wx,  thereby 
yielding  the  acoustic  motion  equation 


d2p 

dx2 


1  d2p 
c2  dt2  ’ 


where  c2 


1 

Ksp' 


The  one-dimensional  wave  equation  above  can  be  easily  extended  to  the  three-dimensional  case. 
Equation  (3.6),  relating  pressure  gradient  to  fluid  acceleration,  becomes 


(P  +  S) 


dw 

~dt 


— Vp 


where  w,  velocity,  is  a  vector.  Equation  (3.7),  relating  velocity  and  compression  of  the  fluid,  is 
now 

d8 
dt 

4[22],  pg.  242. 


P  (1  +  KsP)  V  •  W  -  Kspw  ■  (Vp) 
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If  the  same  assumptions  from  the  one-dimensional  case  regarding  the  significance  of  small  quanti¬ 
ties  are  applied,  the  equations  are  now 


<9w 

“Vp  =  p~m 

85 

Ki-  =  -v-w 

By  taking  the  divergence  of  the  first  equation,  and  solving  the  second  to  eliminate  w  gives 

_  2  d2p  1  d2p 

v-vp  =vP  =  K,pw  =  ^w 


commonly  known  as  the  three-dimensional  wave  equation. 


3.2  Application  to  the  Two  Layer  Waveguide  Problem 


The  approach  to  the  two  layer  waveguide  problem  begins  with  the  partial  differential  equation 
governing  the  propagation  of  waves  through  a  medium,  also  known  as  "the  wave  equation"  in 
homogeneous  form 


V2p- 


1  d2p 
c2  8 12 


0. 


(3.8) 


Time-harmonic  motion  refers  to  the  repetition  of  function  behavior  within  a  specified  period.  Let 
T  (t)  =  Ae~luJt,  where  i  =  y/—l,  so  that  T  is  time-harmonic.  Solutions  to  the  wave  equation  are 
of  the  form  p  =  P  (r,  9,  z)  T  (t).  Then 


V2  (PT) 


1  82  (PT) 
c2  8t 2 


0. 


Recognizing  that  V  does  not  contain  a  t  component  and  the  time  derivative  does  not  concern  P, 


TN2P 


Assuming  time -harmonic  motion, 


i  p&T 

e2  a«2 


0. 


-  3r  (-^2)  P  =  0- 

Simplification  gives 

V2P+^P  =  0 

and  the  substitution  of  the  wave  number  notation,  k  —  results  in 

V2P  +  k2P  =  0.  (3.9) 


Equation  (3.9)  is  known  as  the  Helmholtz  Equation,  and  describes  acoustic  sound  waves  in  tenns 
of  their  pressure,  P.  For  the  remainder  of  this  paper,  no  distinction  will  be  made  between  lower 
case  p  and  upper  case  P,  that  is  P  =  p,  referred  to  as  "pressure". 
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3.3  Fluid  Flow  Through  a  Porous  Media 

According  to  Morse  and  Ingard,5  the  equation  of  continuity  (equation  (3.5))  becomes 

-v  ■  w  =  KpHWt 

for  fluid  flow  through  a  rigid,  porous  solid.  kp  replaces  ks,  and  is  the  effective  compressibility  of  the 
fluid  in  the  pores.  According  to  Morse  and  Ingard,  kp  =  ks  at  high  frequencies,  since  the  behavior 
of  the  fluid  in  the  pores  is  best  modeled  adiabatically.  At  low  frequencies,  the  presence  of  the  solid 
material  tends  to  hold  the  temperature  constant,  thus,  np  =  kt,  the  isothermal  compressibility  of 
the  fluid.6  Conservation  of  momentum  is  also  changed  within  the  porous  medium,  and  becomes 


iupp 


1  Ti¬ 


er 


upv 


w  =  S/p, 


where  a  is  the  flow  resistivity,  and  pp  is  the  effective  density  of  the  fluid  in  the  pores.  pp  is 
expected  to  be  1.5  to  5  times  greater  than  p,  the  original  fluid  density,  because  of  the  tortuosity  of 
the  medium.  The  tortuosity  parameter  measures  the  deviation  of  the  pore  orientation  of  the  medium 

from  uniform.7  A  commonly  used  mathematical  definition  of  tortuosity  is  ,  where  rrid  is 

the  length  of  the  mean  flow  path  between  two  points  and  Ei  is  the  linear  distance  between  the  two 
points.  Let  w  =  —  V'T.  Thenp  =  pp  (d^ /dt)  +  a\l/,  and 


V2\f  =  kpppH- 


n,dV 

+  KvcrEt  —  . 


Dt2 

With  the  exception  of  the  last  term,  this  is  the  wave  equation  derived  in  equation  (3.9),  with  a  wave 


velocity  of  cp  =  J 1/ npppH.  Assuming  time-harmonic  motion,  as  in  the  derivation  above,  gives 


Vz^+  I  —  )  'T  =  0, 

ce 


(3.10) 


where 


l 


Note  that  since  p  =  (—iujpp  +  <r)T'  in  the  time-harmonic  case,  p  also  satisfies  equation  (3.10). 


5 [22],  pg.  253. 

6[22],  pg.  253. 

7[22],  pg.  253. 


Chapter  4 

The  Two  Layer  Waveguide 


The  two  layer  waveguide  problem  is  the  simplest  model  of  the  soil-air  interface  (Figure  4.1).  A 
layer  of  soil  of  height  is  bounded  by  a  rigid  surface  on  the  bottom,  and  by  the  atmosphere  on 
the  top.  The  partial  differential  equation  governing  wave  propagation  in  the  soil  and  air  is  the 
Helmholtz  Equation,  eqn  (3.9),  which  will  be  solved  in  cylindrical  coordinates. 


4.1  Analytical  Approach  and  Separation  of  Variables 

First  consider  the  homogeneous  problem.  Wave  propagation  through  the  air  and  soil,  respectively, 
is  described  through  the  following  two  equations 

V2pa  +  k2apa  =  0 

V2ps  +  k2sps  =  0, 

A2+iwe 

where  ka  =  —  and  ks  =  \  — from  substitutions  into  equations  (3.8)  and  (3.10).  Since  the 

ca  y  Cs 

partial  differential  equations  are  identical  with  the  exception  of  the  subscripts  that  determine  the 
appropriate  boundary  conditions,  a  "generalized  form"  of  the  Helmholtz  Equation  is  used  until 
boundary  conditions  with  respect  to  depth  are  applied.  The  generalized  fonn  of  the  Helmholtz 
Equation  neglects  the  difference  in  acoustic  properties  between  the  atmospheric  and  soil  layers, 
which  do  not  effect  the  general  separation  of  variables  solution  with  regard  to  the  radial  or  angular 
conditions. 

V2p  +  k2p  =  0 

Cylindrical  coordinates  were  selected  due  to  the  cylindrical  features  of  the  membrane  problem. 1 
The  Helmholtz  Equation,  rewritten  in  cylindrical  coordinates,  is 

(r9^)  +  +  d-l  +  k 2v  =  0 

r  dr  y  dr  )  r2  d62  dz 2 

'[16],  pg.  2051. 
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Figure  4.1:  A  mathematical  schematic  of  the  two  layer  waveguide  in  cylindrical  coordinates  dis¬ 
playing  this  paper’s  notation.  pa  and  ps  are  the  density  of  air  and  soil,  respectively,  and  oj  is  the 
angular  frequency  imparted  from  the  point  source. 


Separation  of  variables,  assuming  p  (r,  9,z)  —  F  (r)  G  ( 9 )  Z  (z),  yields 

2_Gy«)  =  2  _  Z^£) 

r  F  (r)  r2  G  (9)  '  Z  (z)  ' 


(4.1) 


Equation  (4.1)  indicates  that  the  two  arguments  on  either  side  of  the  equation  must  equal  the  same 
constant.  This  assumption  is  valid  because  differentiation  of  both  sides  with  respect  to  z  implies 
the  derivative  of  zz^  is  0.  Therefore,  zz^  is  a  constant.  This  argument  works  for  each  variable 
in  equation  (4. 1).  Since  sound  waves  propagate  in  a  way  that  resembles  sinusoidal  functions  (and 
not  exponentials),  choose  —  p  be  an  arbitrary  constant.  Then 


Z"  (z) 
Z(z) 


-k2 


-T 


(4.2) 


and 


1  j  (rf*  M)  1  G"  W 

r  F  (r)  '  r2  G  (9) 


(4.3) 


Continuing  with  separation  of  variables,  consider  equation  (4.3). 
ranging  yields 


fr(rF’(r))  , 


F 


+  r  p 


G"  ( 9 ) 
~GW 


Multiplication  by  r2  and  rear- 


(4.4) 
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4.1.1  Theta  Equation 


Equation  (4.4)  indicates  that  the  two  arguments  must  equal  constants,  from  the  same  argument 
given  for  equation  (4.1).  Choose  A  to  create  sinusoidal  functions  that  will  satisfy  the  periodicity 
conditions  in  equation  (4.8). 


G"  ( 9 ) 


r 


JL 

dr 


(■ rF '  (r)) 
F  (r) 


+  yur2  =  A 


(4.5) 

(4.6) 


Equation  (4.5)  is  equal  to 


G"  {9)  +  AG  ( 9 )  =  0 


(4.7) 


Assume  A  >  0.  (A  <  0  result  in  trivial  solutions)  This  assumption  is  made  based  on  the  periodicity 
condition  implied  by  the  problem  configuration.2  Therefore,  solutions  are  of  the  form 


G  (9)  —  A  cos  ^\/A 9^  +  B  sin 

Boundary  conditions  are  now  imposed,  as  described  by  the  periodicity  condition  below.  The 
solution  G  must  therefore  satisfy  the  periodicity  conditions 


G  ( — 7r)  =  G(tt) 


(4.8) 


G'  ( — 7r)  =  G'  (7 r) 

for  continuity  at  n.  After  imposing  the  boundary  conditions,  the  solution  becomes 

G  (9)  =  A  cos  ( m9 ) ,  (4.9) 


with  A  =  m2  where  m  =  0, 1,  2, .. 


Orthogonality  of  the  Theta  Eigenfunctions 

Orthogonality  is  a  particularly  useful  tool  in  solving  partial  differential  equations.  Because  orthog¬ 
onality  will  be  used  in  later  portions  of  the  mathematical  analysis  of  this  problem,  it  is  necessary 
to  prove  the  orthogonality  of  the  theta-dependent  function,  G.  By  definition,  a  set  of  functions, 
{4>n}™=i  is  orthogonal  on  [0,  L\  if  fQL  4>n  ( x )  4>m  ( x )  dx  =  0,  when  n  ^  m.  Using  the  trigonometric 
identity  cos  (a)  cos  ( b )  =  4  [cos  (a  +  b)  +  cos  (a  —  6)],  we  obtain  that 


A  A 


An  cos  (n9)  Am  cos  (m9)  d9 

r 

[cos  (( n  +  m)  9)  +  cos  (( n  —  m)  9)]  d9 


fill], pg.  307. 
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Performing  the  integral  results  in 

c  _( 

Onm  S 

The  Kronecker  delta  function  results  from  An 
An  forms  an  orthogonal,  nonnalized  basis. 

4.1.2  Radial  Condition 

The  radial  condition  on  pressure  is  established  through  solving  for  F  ( r )  (eqn.(4.6)) 

c) 

r—  ( rF '  (r))  =  —  nr2F  (r)  +  m2F  (r) . 
or 

Rearrangement  and  multiplication  by  r2  results  in 

r2F"  ( r )  +  rF'  (r)  +  (/rr2  —  m2)  F  (r)  =  0.  (4.10) 

Let  (2  =  /rr2.  Then  equation  (4.10)  becomes  a  Bessel  Differential  Equation3 

C 2F"  (0  +  (F'  (0  +  (C2  -  rn2)  F  (0  =  o. 

Solutions  to  the  Bessel  Differential  Equation  are  known  as  Bessel  Functions  of  the  first  and  second 
kind.  Mathematically,  these  are  represented  by  Jm  (C)and  Ym  ((') .  By  the  principle  of  super¬ 
position,  any  linear  combination  of  Bessel  functions  is  also  a  solution  to  the  Bessel  Differential 
Equation.  One  of  these  combinations  is  known  as  a  Hankel  Function  (represented  by  Hm  ((  )), 
which  is  linear  combinations  of  Bessel  Functions  of  the  first  and  second  kind 

h  (C)  =  Hi11  (0  =  Jrn( 0  +  8'.  (0  (4.11) 

h  (C)  =  H‘i>  (0  =  Jm  (C)  -  iYm  (C)  ■  (4.12) 


1  m  —  n 
0  m  n 


=  \  -  in  the  case  where  m  =  n.  This  selection  of 


4.1.3  Depth  Equation 

We  now  consider  the  z  equation,  equation  (4.2),  written  in  terms  of  the  value  of  z. 

-k2a  =  -/i  0  <  z  <  za 
—  k2  =  —/i  zs  <  z  <  0 

Rearrangement,  multiplication  by  —  Z  (z) ,  and  subtraction  of  ( k 2  —  /i)  results  in 


FM. 

Za(z ) 
Z'Az) 
Zs(z) 


3[11],  pg.  306. 


Z”  (z)  +  (kl  -  n)  Za  (z)  =  0  0  <  z  <  za 
Z's  (z)  +  (k2  —  n)  Zs(z)  =  0  zs  <  z  <0 


(4-13) 
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Suppose  >  1 1 .  Solutions  are  then  of  the  fonn 


Z{z) 


A  cos 
C  cos 


(VW^z) 

(\/ks  -  Pzj 


+  B  sin 
+  D  sin 


\Zkl~  Tz) 

Vks~  Ez) 


0  <  z  <  za 
zs<  z  <0 


As  shown  in  Figure  (4.1),  there  are  several  boundary  conditions  that  define  the  solution  in  the  z 
direction.  First,  it  is  necessary  to  define  an  upper  limit,  za,  for  the  atmosphere  when  using  a  point 
source.  At  and  above  this  limit,  the  atmosphere  is  assumed  to  be  a  vacuum,  that  is,  the  pressure  is 
0.  Mathematically,  this  is  written  as 

p(za)  =  0.  (4.14) 

At  the  soil-rigid  substrate  interface,  it  is  expected  that  some  of  the  energy  is  reflected,  while  some 
energy  is  also  absorbed  by  the  rigid  substrate.  However,  the  pressure  at  the  interface  should  be 
constant.  Therefore, 

(zs)  =  0  (4.15) 

At  the  interface  between  the  atmosphere  and  soil,  (z  =  0),  the  atmospheric  pressure  should  be 
equal  to  the  acoustic  pressure  in  the  soil  (Continuity  of  Pressure). 


Pa  (0)  =  ps  (0) 


(4.16) 


Similarly,  at  the  interface  between  the  atmosphere  and  soil,  the  velocity  of  the  acoustic  wave  should 
be  equal  in  both  the  atmosphere  and  soil.  Acoustic  impedance  relates  the  fluid  velocity  to  the 
acoustic  pressure.4  Mathematically,  acoustic  impedance5  is  the  ratio  between  acoustic  pressure,  p, 
and  fluid  velocity,  represented  by  the  following  equation  in  this  problem 


PaU 


(4.17) 


Note  that  this  is  also  continuity  of  velocity  across  the  air-soil  boundary.  Equation  (4.17)  is  derived 
by  assuming  time-harmonic  oscillations  in  equation  (3.4).  Applying  the  boundary  conditions 
(eqns.  (4.14),  (4.15),  (4.16),  and  (4.17))  results  in 


0  <  z  <  za 
zs  <  z  <  0 

where  Bx  is  an  arbitrary  constant.  From  the  boundary  conditions, 


sin  [  \Jk2a  -  pzaj  cos  [y/k2s  -  pzs) 

r  b  i 

r  o  l 

.  (cos '  (sinzs-v/fc2-//).v/fc2_/i 

D 

0 

L  pau  psu  J 

i 

4[22],  pg.  376,  423. 

5 [22],  pg.  259. 


Z(z)  = 


cos  \yjk%  -  fiZsJ  sin  ( sjk\  -  p  (za  -  z) 
\IW~  PZa)  COS  (y/k%  -  11  (z  ~  Z8) 


CHAPTER  4.  THE  TWO  LAYER  WAVEGUIDE 


23 


To  ensure  that  B  and  D  are  non-zero,  the  eigenvalues,  //,  must  force  the  detenninant  of  the  coeffi¬ 
cient  matrix  equal  to  zero.  The  determinant  of  the  matrix,  set  equal  to  zero,  is 


Ps  -\Aa  _  P  COS  VK  -  P  COS  yj  k2  -  pL 

T Pa  VW~ ~P  sin  za \Jh%-  n  sin  y/k 
=  0. 


(4.18) 


s  -  P 


The  solutions  to  equation  (4.18)  must  be  found  numerically.  A  MATLAB™  program  was  written 
to  solve  for  the  eigenvalues.  In  order  to  program  this  specific  function  into  MATLAB™,  it  was 
necessary  to  perfonn  a  change  of  variables.  This  was  done  so  that  the  eigenvalues  would  be 
evenly  spaced,  and  thus,  could  be  found  numerically  over  a  long  range  of  values.  Therefore,  let 

K  =  Zay/k2  ~  /I.  SO 


p  =  K  -  (  - 


K 


and 


Ps  |  cos  (k)  cos  zs\I  k2s  -  k2a  +  (  —  )  ^ 

Zn  J  i  Zn 


(4.19) 


+Pa  I  sin  (k)  sin  zs\j  k2-k2a+[  — 

\  Zn 


K-K+l-l  =0- 

Zn. 


Using  MATLAB 1M  and  a  program  for  the  secant  method6,  eigenvalues  for  the  depth  condition 
were  detennined. 


Air  and  Soil  Eigenvalues 


Note  that  in  equation  (4.18),  pa  <C  ps.  Therefore,  eigenvalues  occur  when 

cos zay]k2  —  /i  cos zs \J k2  —  p  ~  0. 

If  k  is  taken  to  be  za  y/k2  —  //,  then  two  sequences  of  eigenvalues  exist;  one  sequence  exists  for  the 
air  eigenvalues  and  another  for  the  soil  eigenvalues.  The  air  eigenvalues  (na^n)  are  given  by 

«0)n  «  (2 n  -  1) 

while  the  soil  eigenvalues  (ns,n)  are  given  by 


6[17],pg.  70. 


ZsVks  ~  Pn  =  Zsj k2  —  k2+  ^ 


K, 


(2  n 


1) 


7T 

2’ 
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Figure  4.2:  This  figure  shows  the  soil  and  air  eigenvalue  plots,  respectively,  for  the  two-layer 
waveguide. 


or 


-  k2s  +  kl 


for  n  G  N.  Resonances  occur  when  Ks,n  ceases  to  be  evanescent  and  when  KStH  — >  0.  Denote 
these  as  E-resonances  and  Z-resonances.  An  E-resonance  occurs  when  nSjU  =  kaza .  This  results 
in 


Cp,s 


(2  n 


1) 


7 r 


2  k. 


/  Cp,s 


(2 n  -  1) 
4  IzJ 


Similarly,  a  Z-resonance  occurs  when  K,%n  =  0,  which  gives 


fn 


(2 n  -  1) 


C<7,  C  c 


V~c 


The  location  of  the  E-resonances  and  Z-resonances  are  shown  in  Figure  (4.2).  Note  that  in  this 
project,  the  input  values  that  generated  each  graph  are  listed  in  the  appendix. 
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Characteristic  Function  at  50  Hz 


Figure  4.3:  The  first  six  eigenvalues  for  the  depth  condition.  The  k  values  are  limited  from  0  to 
20  in  this  graphic,  for  demonstration  purposes.  The  zeros  of  the  characteristic  function  determine 
the  eigenvalues. 


Eigenvalue  Plots 

Figure  (4.3)  shows  a  plot  of  the  characteristic  function  (eqn.  4.19)  at  /  =  50  Hz.  The  eigenvalues 
in  this  plot  are  the  E-resonances.  Figure  (4.3)  indicates  the  first  of  the  k  eigenvalues  occurs  at 
approximately  2.2,  a  second  eigenvalue  of  5.3,  and  subsequent  eigenvalues  occurring  at  intervals  of 
approximately  ^ r.  A  MATLAB™  program  based  on  the  secant  method  of  finding  zeros  of  a  function 
was  written,  since  the  secant  method  converges  quickly  when  the  interval  containing  a  zero  is 
known.7  The  eigenvalues  depend  upon  the  frequency.  Figure  (4.4)  shows  that  the  eigenvalues  are 
quite  sensitive  to  frequency.  Note,  however,  that  the  spacing  of  the  values  of  interest  for  k  remains 
at  intervals  of  approximately  7r.  Although  there  are  exceptions  to  this  spacing,  the  secant  method 
has  been  highly  effective  in  finding  the  k  which  result  in  zeros  of  the  eigenfunction. 


Orthogonality  of  the  Depth  Condition 


Using  a  Sturm-Liouville  argument,  which  is  outlined  in  Appendix  (A),  the  eigenfunctions  de¬ 
scribed  in  equation  (4.18)  were  confirmed  to  be  orthogonal.  The  weight  function,  as  defined  in 
Appendix  (A),  is 


f  ps  when  0  <  z  <  za 
\  pa  when  zs  <  z  <  0 


7[17],pg.  70. 
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Characteristic  Function  Values 


Figure  4.4:  Eigenvalues  for  the  depth  condition  at  multiple  frequencies,  demonstrating  the  eigen¬ 
function  dependence  on  frequency. 
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4.2  Insertion  of  a  Point  Source 

Insertion  of  a  point  source,  also  known  as  a  unit  impulse,  in  the  upper  waveguide  layer  leads  to  the 
boundary  value  problem 

V2p  +  k2p  —  —  -<5  (r  —  r0)  S  (6  —  0O)  S  (z  —  z0) .  (4.20) 

r 

The  Dirac  delta  "function",  5  (. x ),  is  a  mathematical  representation  of  a  function  equal  to  oo  at 
x  —  0,  and  is  equal  to  0  for  all  other  values.8  In  other  words, 

S(x-Xi)  =  (  °  (4.21) 

I  oo  j:  =  Sj 

Additionally,  the  delta  function  has  the  property  that  the  integral  of  a  function,  /  (x) ,  multiplied 
by  a  delta  function  <5  (x  —  x0 ),  over  the  interval  (a,  b)  (assuming  that  a  <  x0  <  b) 

Jf  (x)  S(x-x0)  =  f  (x0)  ■  (4.22) 

a 

oo  oo 

Let  p  (r,  9,  z)  —  E  E  F mn  (r)  Gm  ( 9 )  Zn  (z).  Then  V2p  in  cylindrical  coordinates  is  repre- 

m= 0  77-—  1 

sented  by 

OO  OO  1 

V2P  =  EE  Kn  (r)  +  -i^B  (r)  Gm  (0)  ZB  (s) 

777, =0  77-—  1  L  ^ 

TL^ 

~~^Fmn  (r)  G,m  (0)  (*)  +  Fmn  (r)  Gm  (6)  [-  ( k 2  -  /i)  Zn  (z)]  . 

Substitution  into  eqn  (4.20)  yields 

oo  oo  1  TT1  ^ 

E  E  C,  M  +  -Fin  (r)  -  —Fmn  (r)  +  »nFmn  (r)  Gm  (0)  Z„  (z) 

m=0 n=l  L  r  / 

=  --5(r  —  r0)  8(9  ^  90)  S  (z -  z0) . 
r 

Multiplying  by  Gk  ( 9 )  y  (z)  Zt  (z)  m,  integrating  to  use  the  orthogonality  of  the  angular  and  depth 
equations,  and  using  the  properties  of  delta  functions  (eqns.  (4.21)  and  (4.22))  gives 

( (r)  +  \f'm  (r)  -  ’Efu  (r)  +  ^Fu  (r))  }  Gk  (9)2  dSfx  M  Z,  (zf  dz 
\  I  '  J  -7T  Zg 

=  --S  (r  -  r 0)  Gk  (0O)  y  (zo)  Zi  (z0)  . 
r 

8[H],  Pg-  392. 
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Letting 


Gk  (flp)  X  {zq)  Zg  (-gp) 

TT  Za 

jGk(9)2d9jX(z)Zl(z)2dz 

—  7T  Zs 


results  in 

rFkl  (r)  +  Fja  (r)  +  (r»n  ~  ^  Fkl  (r) 

The  equation  above  can  be  rewritten  as 

^  (rFh  (r))  +  (rfxn  -  ^  Fkl  (r)  = 

and  integrating  from  r0  —  £  to  r0  +  £  gives 

(r o  +  e)  Fh  (r0  +  e)  -  (r0  -  e)  F^  (r0  -  e)  + 
=  -Cfc, 

Letting  £  — ■>  0  and  requiring  that  Fki  (r)  be  continuous  at  r  = 


=  -CW  (r  -  r 0) . 


- CHS  ( r  -  r 0) . 


U)+£ 

/ 


r  o—£ 


Fu  ( r )  dr 


r0  gives  the  conditions 


(4.23) 


Fu  {t o+)  =  Fh  {t q— ) 


(4.24) 


F'kl  (r0+)  -  F'kl  (r0— )  = - C„ 

Except  at  r  =  r0,  equation  (4.23)  is  just  Bessel’s  equation.  Thus 

n  t  x  /  ^fci-4  (v/E(r)  >  o  <  r  <  r0 
=  n,<r 

is  bounded  at  the  origin  and  satisfies  the  Sommerfeld  Radiation  Condition.9  Bessel  functions  of 
the  first  kind  satisfy  boundedness  at  the  origin.  The  Sommerfeld  Radiation  Condition  states  that 


lim 

r— >oo 


,  Q 

r  (  — w  —  ikw  )  =  0. 
or 


Type  1  Hankel  functions  satisfy  the  radiation  condition.  Additionally,  the  Hankel  function  of 
the  first  kind  indicates  waves  propagating  outward  from  the  source,  which  represents  the  physical 
reality.  The  conditions  expressed  in  (4.24)  require 


Aki'h  (V/V'o)  -  BklH™  (v7¥"o)  =  0 


9[25],  Ch.  6.32. 
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Solving  for  Ay  and  By  results  in 

Ay  = 


By  — 


Jk  iy/Firo)  +  y/WiHk*  iVJpr  o)  J'k  (\//firo) 
~CyJk  (y/JIiro) 


i^/Firo)  Jk  (v^r0)  +  (v'/V'o)  J'k  iVJJjo) 

According  to  the  Handbook  of  Mathematical  Functions10,  particularly, 

-zCm+1  (z)  +  mCm  (z) 


C'm  (z)  = 


and 


HW 

nk+ 1 


(z)H 


(2) 


—  H (2^ 

nk+ 1 


the  denominator  of  equation  (4.25)  simplifies  to 


4  i 

7 TZ? 


Thus 


Noting 


-\fPiHk]'  {VJpr  o)  4  (V/fi'ro)  +  V/P^  (v^firo)  J'k  (V/firo)  = 


2* 


^CyH^  (yTfi^o)  _  77T 


-4fe/  = 


27 

7rr0 


=  (v^r0) 


^CkiJk{y/Wiro)  in  / 

2~ 


Bn  - r"  . "  "  -  =  CkiJk  {Veto)  . 


27 

7rro 


f  G o  ((9)2  dO  =  27t,  when  G0  (0)  =  1 


nr0 


gives 


/  Gm  (0)2  d6  =  7r,  when  Gm  ( 9 )  =  cos  (m0) ,  sin  (md) 


Ct)n  — 


X  (*o)  Zn  (z0) 


r  — 

YVm.n. 


2nfx  (z)  Zn  (z)2  dz 

Zs 

Gm  ( 00 )  X  (zo)  Zn  (zq) 


nfx(z)  Zn(z)2  dz 


,  Gm  ( 6 )  =  cos  ( md ) ,  sin  ( m0 ) 


10 


[2],  pg.  362,  eqn.  (9.1.27). 


(4.25) 

(4.26) 

(4.27) 
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Thus,  for  r  <  r0 

p(rAz)  =  (^r0)j0(^r)X(z0)Zn(z0)Zn(z)  (4.28) 

n=  1  ^ n 

oo  oo  ^ 

+  E  E  iHm  {VJEr o)  Jm  (VKr)  cos  (m0o)  cos  (m0) 

m=l  n=l 

+h£]  (vT^ro)  Jm  (v&)  sin  (m0o)  sin  (m0))  x  (x0)  Zn  (s0)  Zn  (s) 

°°  2  /-.  \ 

=  E  -  T^T #0  (v'/V’o)  <4  (vKr)  X  (so)  Zn  (so)  (s) 

71=1 

oo  oo  T 

+  E  E  ~^rHm  {VKr o)  Jm  (y/Kr)  cos  (m  (00  -  0)) 

m= 1  n=l  SSI/ji 

XX  (s0)  Zn  (S0)  Zn  (z) 

and  for  r  >  r0 

OO  a  ,  v 

j>M,z)  =  vS^(vS'„)xWz.Wz.W  (429) 

71=1  4<^71 

oo  oo  T 

+  E  E  ~TTrHm  (vKr)  (vKro)  COS  (m  (0O  -  0))  X  (s0)  (s0)  Zn  (x) 

771=1  71=1  ^^71 

Z a  ^ 

where  4>n  =  fx  (s)  Z„  (x)  fix. 

-2-s 

Velocity  in  the  x-direction  is  given  by 

w2  =  — —  (d2p)  (4.30) 

LOp 

= - ’(^o  *  (^Zr)x(zo)Z„(z0)Z'n(z) 

UP  |_n=l  41n 

OO  OO  -7 

+  E  E  (v'/V’o)  (v//^r)  cos  (m  (°0  -  0))  X  (so)  Z„  (x0)  Z'n  (x) 

771=1  71=1  ^^71 

when  r  <  r0,  and 

Wz  = - E  -jxrH 0  (VKr)  J0  (v'/V'o)  X  (so)  Z„  (x0)  Z'n  (x) 

up  |_n=l 

OO  OO  7 

+  E  E  (VJEr)  Jm  (vKro)  COS  (m  (00  -  0)) 

771=1  71=1  ^^71 

XX  (so)  Zn  (x0)  z;  (x)] 

when  r  >  r0.  Note  that  these  results  are  comparable  with  similar  models  found  in  the  literature.11 
11  [9],  [6],  [12],  [23], 
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MATLAB™  User  Interface  for  the  Two 
Layer  Waveguide 


A  user  interface,  featuring  an  experimental  and  laboratory  input,  was  created  from  the  equations 
for  the  Two  Layer  Waveguide.  Specific  instructions  for  the  user  interface  are  included,  as  well  as 
the  theoretical  background  for  the  transition  of  the  soil  properties  to  the  porous  representation. 


5.1  Instructions  for  the  Soil  Resonance  User  Interface 

The  soil  resonance  user  interface  is  designed  to  aid  experimental  physicists  in  predicting  soil  res¬ 
onances  in  field  experiments,  where  equipment  is  limited,  and  in  laboratory  settings,  with  more 
precisely  defined  measurements  and  analysis  tools.  The  field  experiment  case  ("basic  module"), 
which  utilizes  atmospheric  temperature  and  humidity  and  relies  on  more  qualitative  measures  for 
soil  type,  moisture,  and  compactness,  is  designed  for  rough  predictions  based  on  previous  data. 
The  laboratory  case  ("advanced  module")  uses  precise  laboratory  tests  to  detennine  the  density 
and  sound  speed  through  the  air  and  soil. 

5.1.1  General  Instructions  for  the  User  Interface 

1 .  Upon  obtaining  the  User  Interface  files,  save  the  files  and  note  their  location. 

2.  Open  MATLAB™. 

3.  In  the  "Current  Directory"  script  next  to  the  help  icon  (a  yellow  ?),  select  the  file  directory 
where  the  Interface  files  are  stored. 

4.  Type  "gui"  into  the  "Command  Window". 
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5.1.2  Instructions  for  Use  of  the  Basic  Module 

1 .  To  select  the  basic  module,  press  enter  when  presented  with  the  advanced  selection. 

2.  Enter  the  ambient  temperature  in  degrees  Celsius. 

3.  Enter  the  ambient  humidity  in  decimal  form. 

4.  Select  the  appropriate  soil  type  which  most  closely  describes  the  soil  in  the  experiment. 
Select  only  the  numbers  given,  decimal  values  are  not  acceptable. 

5.  Approximate  the  soil  moisture  on  a  scale  of  1  to  5,  with  1  being  very  dry  soil,  to  5  being 
very  moist  soil. 

6.  Approximate  the  packing  of  the  soil  on  a  scale  of  1  to  4,  with  1  being  very  loose  soil,  and  4 
being  very  compact  soil. 

5.1.3  Instructions  for  Use  of  the  Advanced  Module 

1 .  To  select  the  advanced  module,  press  any  number  followed  by  enter  when  presented  with  the 
advanced  selection. 

2.  Enter  the  speed  of  sound  in  the  atmosphere.  (Default  -  343  m/s).  To  select  the  default  value, 
press  enter  without  entering  a  value.  Default  values  are  available  for  all  advanced  menu 
entries.  An  explanation  of  the  chosen  default  value  follows  in  the  "Theoretical  Background" 
section  below. 

3.  Enter  the  speed  of  shear  sound  waves  in  the  soil.  (Default  -  2 10.3 1  m/s) 

4.  Enter  the  speed  of  compressional  sound  waves  in  the  soil.  (Default  -  647  m/s) 

5.  Enter  the  flow  resistivity  of  the  soil1 .  (Default  -  0  ) 

6.  Enter  the  density  of  the  atmosphere.  (Default  -  1.205  ^) 

7.  Enter  the  tolerance  of  the  secant  method  computation.  (Default  -  5e-5) 

5.1.4  General  Instructions  Following  the  Basic  or  Advanced  Module 

1 .  Enter  the  depth  of  the  porous  soil  (above  bedrock)  AS  A  NEGATIVE  VALUE.  For  example, 
if  the  depth  of  the  soil  before  hitting  bedrock  is  0.175  meters,  the  depth  would  be  entered  as 
" -0.175". 

2.  Enter  the  position  of  the  source  as  a  vector.  See  Figure  (5.1)  below  for  further  clarification. 

‘[3],  pg-  175. 
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Figure  5.1:  The  physical  representation  of  the  numerical  entry  quantities  for  source  and  position 
of  the  two-layer  waveguide  user  interface. 


3.  If  interested  in  a  frequency  -  pressure  plot,  enter  the  range  of  frequencies  of  interest  as  a 
vector  (ie,  if  a  frequency  sweep  from  100-200  is  desired,  enter  [100,200],  If  a  depth  - 
pressure  plot  is  of  interest,  press  enter.  Then  enter  the  depth  range  of  interest  (in  the  same 
format  as  the  frequency  sweep),  press  enter  and  enter  the  frequency  of  interest. 

5.2  Theoretical  Background 

The  theoretical  background  necessary  to  write  the  programs  devised  for  the  numerical  implemen¬ 
tation  is  described  in  this  section.  The  background  presented  here  is  applicable  for  the  numerical 
implementation  shown  throughout  this  paper. 

5.2.1  The  Basic  Module 

The  ambient  temperature  is  used  to  approximate  the  density  of  the  air  and  the  sound  speed  through 
the  air.  Density  was  determined  using  data  from  E.  Sengpiel2,  and  sound  speed  was  found  using 
the  formula3 


2[24] 

3  [24] 


Ca)T  =  331.4  +  0.6  x  Tarn\)  (°c). 
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Soil  Type  in  User  Interface 

Soil  Type  in  Oelze  Paper 

Sand 

Plainfield  Soil 

Silt 

Sable  Soil 

Clay 

Drummer  Soil 

Organic  Matter 

Adrian  Soil 

Figure  5.2:  This  table  shows  the  user  interface  soil  labels  with  the  soil  labels  from  the  original 
paper  by  Oelze  et  al4 5 6. 


Humidity  also  affects  sound  speed  and  density.  The  atmospheric  density  was  adjusted  by 

Pa  =  Pq,t(1  +  H)/{1  +  1.609  H),4 

where  pa  T  is  the  density  from  the  temperature  data,  and  H  is  the  humidity  percentage.  The  sound 
speed  through  the  air  is  determined  from 

Ca  =  Ca,T  +  0.6  Hca)T5 

The  soil  type  and  composition  entries  are  compared  to  data  collected  by  Oelze  et  al.6  The  soil 
labels  are  contained  in  table  (5.2).  The  exact  percentage  of  composition  is  also  available  in  the 
Oelze  paper.7 

5.2.2  The  Advanced  Module 

The  default  values  for  density  and  speed  in  each  medium  are  from  data  collected  by  Dr.  D.  Keith 
Wilson  and  Dr.  H.  Cudney  in  a  dry,  desert  -  like  environment.  The  default  a  =  0  was  chosen  for 
simplicity.  Finally,  the  tolerance  for  the  secant  method  is  a  reasonable  value  for  the  data. 


4[26] 

5  [24] 

6[8],pg.  792. 
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Chapter  6 

Numerical  Analysis  of  the  Two  Layer 
Waveguide  Problem 


A  MATLAB™  program  was  devised  to  calculate  the  pressure  and  velocity  described  in  equations 
(4.28)  and  (4.30).  The  program  and  its  supporting  commands  can  be  found  in  the  appendix.  There 
were  many  reasons  behind  writing  the  program  for  the  simple,  two  layer  waveguide  problem.  The 
program  will  provide  the  basic  structure  for  future  programs  involving  more  complicated  problems, 
and  will  allow  the  development  of  solutions  to  anticipated  problems.  Additionally,  the  effects  of 
numerical  selection  for  the  atmospheric  height  (za)  and  rigid  soil  depth  (zs)  need  to  be  investigated 
to  analyze  their  contribution  to  the  pressure  and  velocity.  Assumptions  for  the  variables  for  each 
graphic  are  as  follows  (unless  otherwise  stated):  a  =  0,  cs  =  160,  ca  =  330,  =  —1,  za  =  500. 

The  value  for  a  was  chosen  for  simplicity,  cs,  ca  and  zs  are  appropriate  for  the  environment  of 
interest,  and  za  was  computationally  determined. 


6.1  Determining  Resonances  for  the  Two  Layer  Waveguide  Prob¬ 
lem 

Equation  (4.28)  describes  the  pressure  at  a  point  (r,  6,  z).  Resonant  frequencies,  or  resonances, 
are  the  frequencies  at  which  there  are  large  amplitude  oscillations.  Thus,  resonances  are  indicated 
by  large  deviations  in  pressure  caused  by  the  acoustic  wave.  Resonances  in  Figure  (6.1)  occur  at 
/  ~  136,  228,  319, 412  Hz.  In  context  of  the  two  layer  waveguide  problem,  resonances  indicate  that 
the  soil  surface  oscillates  at  maximum  amplitude  at  these  frequencies.  The  direction  of  the  peak 
(up  or  down)  is  not  significant,  it  is  connected  to  constructive  and  destructive  interference.  Should 
landmine  detection  be  attempted  at  low  soil  depths  (before  hitting  bedrock),  Figure  (6.1)  suggests 
that  resonances  detected  in  the  range  of  /  ~  136, 228, 319, 412  Hz  are  probably  soil  resonances, 
instead  of  mine  resonances.  Figure  (6.2)  compares  the  effect  of  deeper  soil  on  resonant  effects  of 
the  rigid  substrate.  It  is  reasonable  that  resonances  existing  at  a  depth  of  1  meter  would  not  be 
as  apparent  or  appear  at  all  in  deeper  soil,  due  to  greater  attenuation  with  distance.  In  fact,  there 
are  more  resonances  in  deeper  soil,  but  these  resonances  occur  at  broader  frequencies  in  shallow 
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Soil  Surface  Pressure  Profile  in  a  Soil  of  Depth  1m 


Figure  6.1:  Graphical  representation  used  to  find  resonances.  The  resonances  of  the  two-layer 
waveguide  occur  from  the  soil  layer,  and  differ  with  depth. 


Comparison  of  Resonances  at  Different  Depths 


Figure  6.2:  Graphical  representation  of  resonances  demonstrating  the  effect  of  soil  depth. 
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Resonance  Program  Comparison 


Figure  6.3:  Comparison  of  two  independantly  written  programs  for  the  two-layer  waveguide  prob¬ 
lem. 


soil.  The  peaks  at  approximately  /  =  450  are  believed  to  be  resonances  which  happen  to  occur  at 
the  frequencies  used  for  the  calculation.  Note  that  the  resonances  in  the  deeper  soil  are  narrower 
than  the  shallow  depth  (zs)  resonances.  To  ensure  accuracy  in  programming,  both  Mattingly 
and  Buchanan  wrote  independent  codes.  There  is  a  slight  disagreement  in  the  curves,  due  to 
MATLAB™’s  computation  of  the  depth  condition,  which  was  written  in  two  different,  equivalent 
ways  by  the  programmers. 


6.2  Selection  of  a 

According  to  Morse  and  Ingard1,  a  is  the  pressure  drop  required  to  force  a  unit  flow  through  the 
material.  If  the  value  for  a  is  too  high,  the  resonances  will  not  appear,  since  the  flow  resistance  is 
too  high.  Similarly,  if  the  flow  resistance  is  too  low,  the  model  will  not  match  the  physical  reality. 
Figure  (6.4)  shows  the  effects  of  low  and  high  values  of  a. 

'[22],  pg.  253. 
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The  Effects  of  Flow  Resistivity  (cr)  on  Pressure 


Frequency  (Hz) 

Figure  6.4:  This  graph  demonstrates  the  effect  of  a  on  pressure.  Note  that  resonances  do  not 
appear  when  a  is  very  high. 


6.3  Selection  of  za 

The  effects  of  the  numerical  selection  for  the  atmospheric  height  (za)  were  considered.  The  ap¬ 
propriate  value  for  za  will  not  show  significant  deviation  from  much  larger  values,  and  its  behavior 
will  be  consistent  with  the  larger  values.  za  needs  to  be  as  small  as  possible  to  maximize  the 
computational  efficiency  of  the  program.  The  receiver  height  versus  pressure  graph,  Figure  (6.5), 
allows  comparison  of  the  different  values  for  za.  After  investigation,  an  appropriate  value  for  za 
is  500  meters.  The  pressures  recorded  at  500  meters  closely  resemble  the  pressures  recorded  at 
800  meters  over  the  span  of  100  hertz,  as  shown  in  Figure  (6.6).  Lower  values  for  za,  as  show 
in  Figure  (6.5),  show  noticeable  deviation  from  za  =  500  m.  For  example,  the  next  lower  value 
of  za  =  200  m  (in  red)  does  not  closely  resemble  the  behavior  of  za  =  800  m,  particularly  at  the 
resonant  frequency,  /  =  128  Hz. 


6.4  Effects  of 

The  considerations  for  the  depth  of  the  soil  is  entirely  dependent  on  the  operating  environment. 
The  rigid  soil  depth  physically  represents  the  depth  at  which  the  porous  soil  layer  interfaces  with 
a  more  dense  (bedrock-type)  soil  that  reflects  most  of  the  acoustic  energy.  This  soil  depth  is 
also  dependant  on  soil  type.  In  Figure  (6.7),  a  dry  soil  environment,  the  resonant  effects  are  very 
distinguished  when  compared  to  Figure  (6.8),  a  moist  soil  environment.  The  effects  of  moisture 
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Determining  an  Appropriate  Value  for  Atmospheric  Height  (Za) 


Figure  6.5:  Graphical  representation  of  the  effect  of  za  on  the  pressure  placed  on  the  mine.  The 
optimal  value  for  za  will  approximate  za  at  much  larger  values,  while  remaining  computationally 
efficient. 


Determining  an  Appropriate  Value  for  Atmospheric  Height  (Za) 


Figure  6.6:  Graphical  representation  of  the  effect  of  za  on  the  pressure  at  z  =  0.  This  graphic 
shows  the  comparison  for  the  hypothesized  optimal  value  for  za  (500  m),  approximating  za  at  much 
larger  values  (800  m). 
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Effect  of  Soil  Depth  (Zs)  on  Resonance  and  Pressure 


Figure  6.7:  Graphical  representation  of  the  effect  of  zs  on  the  pressure  placed  on  the  mine.  This 
graphic  shows  the  comparison  for  depth  of  the  bedrock  layers  in  a  sandy  soil  environment. 


are  incorporated  into  the  soil  sound  speed,  cs.2  Smaller  pressures  are  recorded  in  the  moist  soil 
environment,  and  the  resonances  are  barely  distinguishable.  This  mathematical  finding  makes 
physical  sense:  a  moist  soil  is  more  dense,  causing  more  attenuation  in  the  sound  wave,  which 
results  in  lower  pressure  differences  and  decreased  resonant  effects,  a  remained  at  zero  for  the 
computation  presented  below,  but  future  analysis  of  different  soils  will  need  to  include  different 
cr,  which  is  clearly  effected  by  soil  type.  Values  suggested  in  a  soil  sound  speed  paper3  are  being 
examined  for  their  relevance  to  this  problem. 


2[8],pg.  792. 
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Effect  of  Soil  Depth  (Zs)  on  Resonance  and  Pressure  (Moist  Soil) 


Frequency  (Hz) 


Figure  6.8:  Graphical  representation  of  the  effect  of  zs  on  the  acoustic  pressure.  This  graphic 
shows  the  comparison  for  depth  of  the  bedrock  layers  in  a  mixed  soil  environment.  (The  compu¬ 
tations  in  this  graphic  considered  only  a  change  in  sound  speed,  since  the  values  for  a  have  yet  to 
be  validated. 


Chapter  7 

Solution  of  the  Membrane  Problem 


Figure  (7.1)  represents  the  membrane  problem,  which  consists  of  the  top  plate  of  the  landmine 
imbedded  in  the  rigid  substrate  from  the  waveguide  of  the  two  layer  waveguide  problem.  While  the 
functions  outside  the  imaginary  cylinder  remain  identical  to  those  from  the  two  layer  waveguide 
problem  (see  eqns.  (4.28)  and  (4.29)),  the  equations  within  the  cylinder  will  be  derived  in  the 
following  2  chapters. 


7.1  Definition  of  the  Green’s  Function 

Suppose  now  that  a  membrane  occupies  the  region  D  =  {(r,  9)  :  r  <  a}.  The  equation  for  the 
motion  of  a  damped  membrane  subject  to  an  external  pressure  p  (r,  0.  t)  is 

TV2u  -  (3dtu  +p  =  pdttu. 

where  T  is  the  tension  of  the  membrane  and  p  is  its  density.  If  the  external  pressure  is  time- 
harmonic  p  (r,  9,t )  =  p  (r,  9)  e~lult  and  the  membrane  is  fixed  at  the  boundary,  then 

V2«  +  k2u  =  — 


u  (a,  9)  =  0 

with  k2  =  ( pu 2  +  idijj)  /T.  Let  Q  (r,  9 ,  r' ,  9’)  denote  the  Green’s  function  for  the  problem,  that  is, 


V\rfi)G  +  k2G  =  —S  (: r  -  r')  6  ( 9  -  9’) 


(7.1) 


G  (a,  9,  r\  9')  =  0,  g  (r,  9,  r’,  9’)  =  g  (/,  9\  r,  9) 


42 


CHAPTER  7.  SOLUTION  OF  THE  MEMBRANE  PROBLEM 


43 


Figure  7.1:  A  mathematical  schematic  of  the  membrane  problem  in  cylindrical  coordinates  dis¬ 
playing  this  paper’s  notation.  T  is  the  tension  of  the  membrane,  which  will  be  examined  for  the 
case  of  the  plastic  landmine. 
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Note  that 

fj[G  (r,  9 ,  r',  9')  V2«  (r,  0)  -  uV\rfi)g  (r,  0,  r',  0')]  rdOdr 

0  — 7T 

=  /  /  -g  (r,  o,  r',  9')  (k2u  (r,  9)  +  (r,  0)^ 

+  ^&2(/  (r,  0,  r',  9')  +  -5  (r  —  r')  5  (9  —  9 ')^  u  (r,  0)  rdOdr 
=  f  f  Q  (r,  9 ,  r1, 9')  p  (r,  9)  rd9dr  +  u  (r1,  9') . 

Since  V  •  (uVv)  =  Vu  ■  Vu  +  uV2v,  the  divergence  theorem  gives 


f  f  [g  (r,  9,  r‘,  9')  V2u  (r,  9)  -  uVfr>0)g  (r,  9,  r',  9’)]  rd9dr 

0  —7V 

=  ff  [v  •  (gvu)  -  vs  •  v«  -  (v  ■  (uvg  -  wg  •  v«))]  dA 

D 

=  /  [gvu  ■  n  —  uTjg  ■  n]  ds  =  0 


in  view  of  the  boundary  conditions.  Thus,  the  solution  to  the  non-homogeneous  membrane  prob¬ 
lem  is 

u  (r,  9)  —  f  g  (r,  9,  r' ,  9')  p  (r',  9')  r'dO'dr'  (7.2) 

'  n  — 


7.2  Determination  of  the  Green’s  Function 

It  remains  to  find  the  Green’s  function.  In  polar  coordinates, 

-< dr  (■ rdrg )  +  ~^deeg  +  k2g  =  0 

except  at  (r,  9)  =  (r',  9').  Seeking  solutions  of  the  fonn  0  (r,  9)  =  F  (r)  G  (9)  gives 

HrF'jr))'  1  G"{9) 

F  (r)  r2  G{9 ) 

r  (rF1  (r))'  ,„2.„2  G"  (9) 

F  (r)  *  G  ( 9 ) 

From  equations  (4.7)  and  (4.8)  the  problem 

G"  ( 9 )  +  AG  (0)  =  0 
G  ( — 7r)  =  G  (tt)  ,  G'  ( — 7r)  =  G'  (tt) 
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has  the  solution 

Am  =  m2,m  =  0,1,2,... 

Go  =  1,  Gm  ( 6 )  =  cos  ( mQ ) ,  sin  (m9) . 

The  problem  for  F  (r)  is 

r  ( rF '  (r))'  +  ( k2r 2  —  m2)  F  ( r )  =  0 
F  (a)  —  0 

Substituting  a  solution  of  the  form  Q  (  r,  6)  =  Y^m=o  Fm  ( r )  Gm  ( 6 )  into  the  Laplacian  gives 

^<9r  (rdrQ)  + 

00  1  ^7  1 

=  £  ;  *  (>•<.  (0)  G,„  m  +  -2Fm  M  <?;  («) 

m= 0 

00  /  1  \  rr,2 
=  X  (  F".  W  +  ('0  G„.  (9)  -  (r)  Gm  (9) 

m=0  '  ' 

Substituting  this  into  (7.1)  gives 

(  Fm  (r)  +  -F'm  (r)  j  Gm  ( 6 )  -  —  Fm  (r)  Gm  (6)  +  k2 Fm  (r)  G,m  (0) 

m= 0  '  ' 

=  -U^r -r')5{6 -6') 


E  (r)  +  F'm  (r)  +  Fm  (r )'j  Gm  (0)  =  -5  (r  -  r')  5  (0  -  0') 

Multiplying  by  Gm/  (0)  and  integrating  gives 

(r)  +  F*m,  (r)  +  (Vr  -  ^  Fm.  (r))  f  Gm,  (0)2  dd 
=  -)(r-r')Gm-(0') 

rF. ",  (r)  +  F'm,  (r)  +  ^fc2r  -  Fm,  (r)  =  (r  -  r') 

where 

r*  (0O 

- • 

./  Gm>  {O'fdO 


Integrating  from  r0  —  £  to  r0  +  £  gives 


(F  +  £)  F'm>  ( r'  +  £)~  ( r'  ~  £ )  F'm>  ( r *  ~£)  + 


rr'+e 


m 


/  2 


Arr - - —  )  Fm>  (r)  dr 


—  Dm' 


(7.3) 
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Letting  s  — >  0  and  requiring  that  Fmt  (r)  be  continuous  at  r  =  r'  gives  the  conditions 


Fm'  (r’+)  =  Fm,  (r'~)  (7.4) 

Hr  (>■'+)  -  F'u,  (r'~)  = 

Except  at  r  =  r',  (7.3)  is  just  Bessel’s  equation.  Since  the  solution  must  be  bounded  at  r  =  0, 

F  (  v  f  An'Jrn'  ( kr )  o  <  r  <  r' 

m  ^  \  (kr)  +  Cm’Ymt  (kr)  r0  <  r  <  a 

The  Hankel  Function  is  replaced  by  its  identity,  the  sum  of  Bessel  Functions  of  the  first  and  second 
kind,  for  simplification  of  adherence  to  the  boundary  conditions.  The  jump  conditions  and  the 
boundary  condition  gives 


An' ’k rr/  (kr  )  (kr  )  CmAn'  (kr  )  —  0 


kAm,J'm,  (kr')  -  kBm,J'm,  (kr')  -  kCmY'm,  (kr')  =  ~Dm, 
Bill1  'kn:  (kci)  4~  CmiYmi  (kci)  0 

Thus,  the  system  to  be  solved  is 


Jm>  (kr') 

-Jm>  (kr') 

-Ym,  (kr')  ' 

Am' 

0 

J’m'  ( kr ') 

-J'm'  (kr') 

-YL'  (kr') 

Bm' 

= 

1 

kr' 

0 

Jm'  (ka) 

P 

A 

Cm' 

0 

The  determinant  of  the  coefficient  matrix  is 


Jm'  (kdj) 


Jm'  (kr')  -Ym<  (kr') 
J'm’ikr')  —Y'n,  (kr') 


T Fm'  (kci) 


Jm'  (kr  )  Jm'  (kr  ) 

JL’(kr')  -J'm,(kr') 

=  -Jm'  (ka)  (-Jm'  (kr')  Y'm,  (kr')  +  J'm,  (kr')  Ym>  (kr')) 


(7.5) 


From  equations  (4.26)  and  (4.27): 


-kr'Ymi+ 1  (kr')  +  m'Ymi  (kr') 


kr' 


Jm’  (kr  ) 


=  Ym>+ 1  (kr')  Jm'  (kr)  - 


+Ym>  (kr') 

m'Ym'  (kr') 
kr' 


-kr' Jm'+i  (kr')  +  m!  Jm>  (kr') 
kr' 


Jm'  (kr')  -  Ym>  (kr')  Jm'+\  (kr')  +  Ym<  (kr')  m  Jm  ^  ^ 


kr' 


—  Ym'- )-i  (kr  )  Jm'  (kr  )  Ym>  (kr  )  Jm'+i  (kr  ) 


CHAPTER  7.  SOLUTION  OF  THE  MEMBRANE  PROBLEM 


47 


By  Abramowitz  and  Stegun1 


Jy+ 1  (%)  Yv  (z)  Ju  (z)  (^) 


7r(^) 


so  our  equation  yields: 


7 r  (At') 


Thus,  by  Cramer’s  Rule 


-4m'  — 


Similarly, 


7r  (At') 


2  Jmi  ( ka ) 


0  -  JTO'  (At')  -Ym/  (At') 

0  Jm>  (ka)  Ymi  (ka) 


*  (hr1)  f _ l_D 

2  Jmi  (ka)  \  kr'  m 


-Jm>  (kr')  -Ym,  (kr') 
Jm'  (ka)  Ym'  (ka) 


7T  (. kr ') 


Noting  that 


2 Jm>  (ka)  \kr ' 


-  .  Dm ' 


■ Jm '  (kr  )  Ymi  (ka)  T-  Ym>  (kr  )  Jm'  (ka)) 


7 r 


2  Jmi  (ka) 


(Dm1)  (  Jm'  (kr  )  Ym/  (ka)  E  Ym /  (kr  )  Jm'  (ka)) 


7r  (At') 


m  2Jm/  (ka) 


7T 


Cm/ 


2Jm/  (fco) 
7r  (At') 


Jm/  (At')  0  (At') 

J'm1  (kr')  -±,Dm,  - Y’m,(kr ') 
0  0  Ym/  (fco) 

D  rn<  J mf  (kr  )  Y^'  (ka) 


2  Jmt  (ka) 


Jm'  (kr')  -Jm'  (kr')  0 

JL'(kr')  -J'm,(kr')  -J.Dm, 


0  Jm'  (ka) 


0 


7 r  7T 

2T  (ka)^>m'Jm'  ^  ^ ^  ^  2  kJm'Jm '  (kr  ) 


Do  = 


Go  (O') 


Dm 


f  G0  (O')2  dO 

— 7T 

Gm  (O')  Gm  (O') 


1 

27T 


/Gm  (0')2^ 


7T 


,  Gm  (0')  =  cos  (O') ,  sin  (0') 


1  [2],  pg.  362,  eqn.  (9.1.16). 
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gives 


0M,r',0')  =  F0(r)  +  J2Fm(r)Gm(d) 


m=  1 


4  J0  ( ka ) 


E 


(J0  (At')  Y0  ( ka )  -  J0  (/ca)  F0  (At'))  Jo  (kr) 

1 

( Jm  ( kv  )  ym  (ka)  Jm  (ka)  Ym  (kv  )) 


^  (  2  J,„  (ka) 
x  Jm  (At)  cos  (m$)  cos  (m#7) 


( Jm  (Arr7)  Em  (fca)  -  Jm  (fca)  (A:/)) 


^  2  Jm  (ka) 

m=  1  v  7 

x  Jm  (kr)  sin  (m$)  sin  (m07 ) 

,  jl(l  t  (Jo  (A:r')  Eo  (Asa)  -  J0  (A;a)  F0  (At'))  J0  (At) 
4  J0  (rea) 


E 


^  2  Jm  (ka) 
x  Jm  (At)  cos  (m  (0  —  (9')) 


(Jm  (At  )Ym  (ka)  Jm  (ka)  Ym  (kr  )) 


for  r  <  r'  and 


0(r,e,r',6')  =  F0(r)  +  J2Fm(r)Gm(6) 


m=  1 


,  J  !,  T  Jo  (At')  *o  (M  Jo  (At)  -  ^  J0  (Air7)  E0  (kr) 
4  J0  (reo)  4 


+  ^  2  Jm  (A:a) 

x  Jm  (kr')  cos  (m  (0  —  07)) 


(Eii  (ka)  Jm  (kr )  Ym  (kr )  Jm  ( ka )) 


for  r  >  r' . 
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Returning  to  (7.2) 


u(v,9 )  =  ^=f  f  G  (r,  9 ,  r',  9')  p  (r7,  9')  r'dO'dr 7 

0  — 7T 

=  f  A  T  7bn)J°  ^  F°  ^  “  Z  J°  ^ 

t  0-5T  L4J0  (fto)  4 

X  1 

T  y  yz  r  (Ym  (ka)  Jm  (kv )  Ym  (kv )  Jm  ( ka ))  Jm  (kv  ) 

i  2  (fca) 

m=l 

x  cos  (m  (0  —  07))]  p  (r7,  07)  r'dO'dr' 

+  bj  I  T  T  lb  \  ( J°  ^r/)  y°  (fca)  _  J°  (fc°)  y°  (^r/))  J0  (Ax) 

t  0-5T  L4^0  (fcQ-) 

X  ^ 

T  ^  y~  r  (Tm  ( kv  )  Zm  (ka)  Jm  (ka)  (kv  ))  Jm  (kv) 

“  2  Jm  (ka) 

m= 1  x  7 

w  Z™  /d  tf'Wl  ^  rt  Art  Art 


x  cos  (m  (0  —  07))]  p  (r7,  07)  r'd9  dr‘ 


7.3  Comparison  with  Mathews  and  Walker 

Mathews  and  Walker2  approach  the  similar  non-homogeneous  membrane  problem 


Vfr,0)G  +  k2G  =  -5(r-r’)5(9-9') 


with  jump  boundary  conditions 


Am’  Jm'  (kv  )  Bm'Jm'  (k'V  )  Cm’Ym'  (kv  )  —  0 


kAm,J'm,  (kv')  -  kBm,J'm,  (kv')  -  kCm.Y^,  (kv')  =  -Dri 


The  Green’s  Function  is  then 


Bm'  J m1  (ka)  T  Cm>Ymf  ( ka )  —  0. 


J]  AmJm(kv)  cos  (m0)  (v  <  r') 

m 

Y  Bm  [ Jm(kv)Ym(ka )  -  Ym  (kv)  Jm  (ka)}  cos  (m9)  (v  >  v') 


where 


Am  =  - -  ,  r  [Jm  (ka)  Ym  (kv')  -  Ym  (ka)  Jm  (kv')} 

^mAm  (Ka) 


B  in. 


Jm  (k’V  ) 

2cm  J rn  (ka) 


2 [20],  pg.  273. 
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and 


Therefore,  for  r  <  r'. 


f  2  if  m!  =  0 
\  1  if  m!  >  0 


Q  =  Yh  AmJm{kr )  cos  (md) 

m 

=  4  t  IT  \  [ J°  ( ka )  Yo  {hr')  -  Y0  (. ka )  J0  (/cr')]  J0(^) 

4c/Q  yrCCL J 

+  E  7  \  [<fm  (&a)  ym  (/cr')  -  (/ca)  Jm  (/cr')] 

m  m  \K>Cl) 

xJm(kr )  cos  (md) 


while  for  r  >  r', 

Q 


E  [Jm(/cr)ym(/ca)  -  Fm  (/cr)  Jm  (ka)]  cos  (m0) 

m 


Y [</0(^)yoM  -  E  (Ax)  J0  Ml 

4t/Q  yrCd J 

+  E  r?  [Jm(/cr)Fm(/ca)  -  Fm  (/cr)  Jm  (/ca)]  cos  (m0) 

m  ^ 'J m  \Kd J 


The  solution  in  Mathews  and  Walker  is  identical  to  the  solution  presented  in  this  paper,  with  the 
exception  of  a  negative  sign  in  each  portion  (due  to  the  absence  of  the  negative  sign  preceding  ^  in 
equation  (7.7)). 


Chapter  8 

Membrane  Imbedded  in  a  Rigid  Substrate 


Let  D  =  {(r,  9)  :  r  <  a}  be  the  area  occupied  by  a  membrane  with  a  point  source  located  at 
(r0,  0o,  z0),  a  <  r0  and  an  atmosphere  satisfying  a  pressure  release  condition  at  z  =  za.  From 
equations  (4.28)  and  (4.29),  the  incident  pressure  is 


°o  oo  ,  , 

E  E  (vKro)  y^r) 

m= 0 n=  1 

x  cos  (m  (6*0  -  6*))  X  {zq)  Zn  (z0)  Zn  (z) , 

OO  OO 

E  E  -sk®.1  (VivO  (vsk 

m=0  n=l 

„  x  cos  (m  (6*0  -  6*))  y  (,s0)  Zn  (z0)  Zn  (z) , 


where 


f  2  m  —  0 

(1  m  >  0 


Total  pressure  is  the  sum  of  the  incident  pressure  field,  pinc,  and  the  scattered  pressure  field,  pscat, 
thus  p  =  pinc  +  Pscat-  Pscat  must  also  satisfy  the  boundary  value  problem 


V2p  +  k2p  =  0 
d 

—p  (r,  6*,  zs)  =  0,  r  >  a 
oz 

p  ( r ,  9 ,  za)  =  0. 
p  (r,  9,  0+)  =  p  (r,  9,  0-) 

,  §~ZP  {r,  9 ,  0+)  =  ,fzp  {r,9,  0-) 

Pau  Ps^ 

Again,  equation  (8.4)  is  derived  from  the  time-harmonic  assumption  on  the  conservation  of  mo¬ 
mentum  equation  (eqn.  (3.4)).  Solutions  are  of  the  form 


(8.1) 

(8.2) 

(8.3) 

(8.4) 


Pscat  (r,  9,  z)  =  £  £  (y/fjQr)  (■ Cmn  cos  (■ m9 )  +  Dmn  sin  (m6*))  Zn  (z) 

m= 0 n— 1 

r  >  a, 
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where  fin  are  the  eigenvalues  from  Problem  two  layer  waveguide  (eqn.  4.18).  Above  the  mem¬ 
brane,  total  pressure  is  given  by 

OO  OO  /  _  \ 

P  (r,  9,z)  =  EE  Jm  (  )  (Ann  cos  (m9)  +  Bmn  sin  (: m6 ))  cpmn  (z) ,  (8.5) 

m= 0  n=  1  '  ' 

r  <  a, 


where  the  eigenvalues  (mn  and  the  eigenfunctions  (pmn  (z)  are  to  be  determined. 

At  the  interface  r  =  a  the  continuity  conditions 

Pscat  (a+,  9,  z)  +  Pine  (a+,  9,z)  =p  (a-,  6,  z)  (8.6) 

drPscat  (a+,  9,  z)  +  d rpinc  (o+,  9,  z)  =  d rp  (a-,  9,  z) 

are  imposed.  At  z  =  zs  continuity  of  velocity  gives,  in  the  case  of  the  time-hannonic  oscillations, 

— iuu  (r,  9)  = - dzp  (r,  9,  zs) ,  r  <  a  (8.7) 

iups 

where  u  is  the  vertical  displacement  of  the  membrane  and  ps  is  the  density  of  the  soil. 


8.1  Solution  of  the  Problem 

8.1.1  The  Eigenvalue  Problem 

The  eigenfunctions  in  the  atmosphere  (i)rnn  a  (z)  satisfy  the  equations 

</>mn,a  (Z)  +  ( kl  ~  C)  Kn,a  (Z)  =  0 

0 mn,a  0 

The  differential  equation  has  solution 

4>mn,a  {z)  =  E  COS  (fiaz)  +  F  Sin  (Oaz)  , 


where 

^  ~  C- 

The  pressure  release  condition  (eqn.  (8.2))  gives 

E  cos  {Daza)  +  Fsin  (Oaza)  =  0 
E  cos  (: Qaza )  =  -F  sin  (: Qaza ) . 

Note  that  (j)mna  (z)  becomes  E1  sin  (Oa  (z  —  za))  for  some  E{.  The  eigenfunctions  in  the  soil 
satisfy  the  equations 

(Z)  +  {ks  -  C)  <l>mn,s  (-)  =  0 
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The  differential  equation  has  solution 

< Kn,s  (z)  =  K  cos  (Msz)  +  L  sin  (Dsz) , 


where 

=  k2s  -  c 

Note  that  equation  (8.1)  does  not  apply,  since  for  <pmn.s,  r  <  a. The  atmosphere-soil  interface 
condition  (eqn.  (8.3))  on  pressure  implies 

-Ei  sin  (Oa  (za))  =  K 


and  continuity  of  velocity  (eqn.  (8.4))  gives 


. Oa  (cos  (Oa^a))  ,DSL 

-i - =  —  i - 


PaN 


PsU 


Substitution  into  (j)mn  a  (z)  yields 


and  (f)mn„  (z)  is 


(*)  =  Ei 


(Pmn,a  (Z)  =  El  sin  (fi0  (z  ~  Za)) 

^ s  a  cos  ( Daza )  sin  (Qsz)  —  sin  (f2az0)  cos 


'mn,s  \~  /  “-l  \  q 

Thus,  the  eigenfunctions  have  the  fonn 

f  sin  (Oa  (z  -  za))  0  <  z  <  za 

<Pmn  (z)  |  cos  (£iaZaj  sin  (OsZ)  -  sin  (Oa  (za))  cos  (Qsz)  zs  <  z  <  0 

From  (7.6)  the  response  of  the  membrane  to  an  external  pressure,  pe,  is 


u(r,9)  =  V -  1  (Ym  ( ka )  Jm  ( kr )  -  Ym  (. kr )  Jrn  ( ka ))  (8.8) 

0  — 7r  m=0  £mJm  \Ka) 

x  Jm  ( kr ')  cos  (m  ( 9  —  9'))pe  (r' ,  9')  r'dB'dr' 

1  an  °°  1 

+  I  c  T  7hn\  ( Jm  ^  Ym  ^  ~  Jm  ^  Ym  ( kr' )) 

r  —TV  n  ^TTl'J HI  \lvQ') 

m— 0 

x  Jm  (kr)  cos  (m  (9  —  9'))  pe  (/,  9')  r'dB’dr' 

1  °o  1 

=  AT  ^  >  t  77  7  (kn)  Jm  (kr )  Ym  (kr )  Jm  (ka)) 

2 T  '  £mJm  (ka) 

m= 0  x  7 

x  J  Jm  (kr')  f  cos  (m  (9  —  9'))  pe  (r1 ,  9')  dO'r'dr ' 

0  —TV 

+  ~T nr \ /  (■ Jm  {kr')  Ym  (ka)  -  Jm  (ka)  Ym  (kr')) 

n  (rvOj)  r 

m— 0  x  7 

x  f  cos  (m  (9  —  9'))  pe  (r1,  9')  dB’r'dr'. 


—  TV 


CHAPTER  8.  MEMBRANE  IMBEDDED  IN  A  RIGID  SUBSTRATE 


54 


Since  in  this  problem  pe  (r',  O')  =  p  (r',  O',  zs )  where  p  is  given  by  (8.5), 


f  cos  (m  (0  —  0'))  pe  ( r O')  dO ' 


7T  OO  OO  /  _  \ 

=  J  Jm>  ( Vc^nr' )  {A-m'n  cos  (m! O')  +  Bmfn  sin  (■ m! 0 ')) 

— 7T  m'= 0  n=l  '  ' 


X(pmn  (zs)  cos  (m  ( 0  —  0'))  dO' 


OO  /  _  \ 

=  £mTT  Yj  Jm  (  )  ( Amn  cos  (5 mO )  +  Bmn  sin  (mO)) 

n=  1  V  / 


x  cos  (: Qaza )  cos  (fiszs) - sin  (Daza)  sin  (fis£s) 

V  Pops  ) 

in  which  case  equation  (8.8)  becomes 

OO  OO  ^ 

u(r,0)  =  -VV  (Ym  ( ka )  Jm  (kr)  -  Ym  ( kr )  Jm  (ha)) 

Z1  m=  0  n=l  ^ a) 


x  (Amn  cos  ( mO )  +  Bmn  sin  (mO))  f  Jm  (kr')  Jm  (^y/(r r'dr' 


cos  (fia^a)  sin  (Oszs)  -  sin  (Daza)  cos  (nszs 


7T  AL  AL  T  (kr) 

+  9T  5]  51  FTT  (  (Amn  cos  (m0)  +  Bmn  sin  (m0)) 

m=0  n= 1  W 

x  f  (Jm  (kr')  Ym  (ka)  -  Jm  (ka)  Ym  (kr'))  Jm  (y/C^r')  r' dr 


cos  (fiaza)  sin  (fiszs)  -  sin  (Ctaza)  cos  (fiszs 


According  to  MAPLE1 


(Yrn  (ka)  Jrn  (kr )  Lm  (kr)  Jm  (ka))  J  Jrn  (kr  )  Jt 


m  \  V  'smn 


r' )  r'  dr' 


+Jm  (kr)  f  (Jm  (kr')  Ym  (ka)  -  Jm  (ka)  Ym  (kr'))  Jm  (  y/C^r'j  r' dr' 

k  f  r  Jm.  (ka)  Jm  (^CYr)  [Jm  (kr)  Ym+1  (kr)  Jm-\- 1  (kr)  Ym  (kr)} 
k2  Cmn  \  EaJm  (kr)  Jm  ( \J Qmna)  [t/m+i  (ka)  Ym  (ka)  Jm  (ka)  Ym_f_i  (Zen)] 

— ~Tj~2  _  >■  7  (^Jm  (kr)  Jm  ^\ZCmna^  ~  ^ rn  (ka)  Jm  ^\/C  mnT^j  ^ 


*  (k2  ~  Cmn) 

where  the  identity1 


‘[2],  pg.  362,  eqn.  (9.1.16). 


Jm+l  (z)  Ym  (z)  Jm  (z)  Ym+i  (z)  — 

7 TZ 
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has  been  used.  Thus, 


^  OO  OO  1 

u(r,6)  =  -V  V  ~fT~2 — 7 — \  t  71  ^  ( Amn  cos  (m0)  +  Bmn  sin  (m0) ) 


X  ^  Jm  {kr )  Jm  ^  \/ C mn^j  Jrn  {kci)  Jm  ^ 

w  ( Ps^a _ ,,  /0  \  „•  /0  \  „•  /o 


^cost^Jstatn.z.l-sin^Jcos^,,) 


and  from  equation  (8.5) 


OO  OO  /  _  \ 

c>2p  (r,  6*,  zs)  =  ]T  ]C  Jm  (  \ZCw~)  (Jmn  cos  (mO)  +  Bmn  sin  (mO)) 

m=0  n=l  '  ' 


x  (  ^  cos  (fiaza)  COS  (0SN)  +  Qs  sin  ( Qaza )  sin  (Q.s2i.s 

V  Pa 


and  equation  (8.7)  becomes 


2  00  00  ^ 

E  E  (p  _  (  )  ,J  (ka)  ('4””  C0S  {m0)  +  S”“”  Sin  (m(,)) 

m=0  n=  1  v  Smn'  m  V  ' 

x  (jm  (hr) 

mn Jm  {kaPj  Jm 

rrivT^  ^ 


=  EEA 


ra=0  n=l 


a  cos  (aA) sin  (sv*)- sin  (sv“) cos  (az-) 


r)  (slmn  cos  (mO)  +  Bmn sin  (m#)) 


m  \  V  Smn' 


x  (  — — -  cos  (f2a£a)  cos  (fE-^)  +  Ds  sin  ( Daza )  sin  (fE^s) 
V  Pa 


Equating  the  coefficients  of  the  trigonometric  functions  gives 


PaN  \  ^  _ 1 _ 

T  “  {k2  -  (mn)  Jm  (ka) 


xAmn  ^  J/n  (Er)  Jm  ^  y/ C mn^j  J<n  ( k(l )  Jm  ^ \/ CrrmT 


x  (  cos  (fia£a)  sin  -  sin  (Daza)  cos  (f^s) 


S  Jm  C mvT  J  Amn 

x  (  cos  (flaza)  cos  ( Qszs )  +  Qs  sin  ( Qaza )  sin  (fEzs 
V  Pa 
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with  the  same  equation  for  Bmn.  Multiplying  by  rJm  ( \J (mnr)  and  integrating  from  0  to  a  gives 


>*U  V  _ 1 _  A 

T  hiV-CmJJmika)  mn 

<  1  l-JLFT  cos  sin  (Dszs)  -  sin  (Oaza)  cos  (Oszs) 

\  oJ-h 


x  .7 


aj  Jr  Jm  ( kr )  J„ 


/’)  dr  —  Jm  (ka)  Jr  J2m  (x/C^r)  dr 


OO  CL  ,  x 

=  E  ( \KZiA  dr  A 

n= 1 0  V  ' 


x  (  cos  (Oaza)  cos  (0.sy5)  +  Os  sin  (Oaza)  sin  (0.sys) 

V  Pa 

Requiring  the  Amn  to  be  arbitrary  gives  the  characteristic  equations 

PsU2 _ 1 _ 

~  T  (k2  -  (mn)  Jm  (ka) 

x  (  COS  (^aZa)  sin  sZs )  ~  sin  (DaZa)  COS  (OsZs)  j 

\Pa*  's  J 


X  J, 


CL  /  _  \ 

=  5rJm  (VC^nr)dr 

n  '  ' 


^a)  }rJm  ( kr )  Jm  (VcEr)  dr  -  Jm  (ka)  JrJ2n  dr) 


:  cos  (Oa)  cos  (SV.)  +  0,  sin  (Oa)  sin  («.*.))  , 


a  _  a 

for  m  =  0, 1,2,3...  Define  h  =  JrJm  (kr)  Jm  (v/CEr)  dr  and  I2  =  frJ2n  (v/CEr)  dr  = 

o  o 

T  \Jrn  +  Jm+ 1  (VCE«)]  •  According  to  MAPLE™, 


h  =  JrJm  (kr)  Jm  (^C~r) 

7  To  V  C mndm  (ka)  Jm-\- 1 

Cmn ~  k2  L 


(/ca)  ^ \/ C mnfl^j  kJm+ 1  (A;a)  Jm  ^\/ 


/2  =  ]V'E 


m  IV  Smn' 


r )  dr  =  —  J„ 


"I”  J m-\- 1 


^m  (  V  Cmn^ 
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which  results  in 


Pau 


T  {k2  -  Cmn)  Jm  (ka) 


X 


Pa^s 


cos  ( Daza )  sin  (Qszs)  -  sin  ( Daza )  cos  (Dszs 


X  ^  Jm  ^\/C iim®|  ^  1  'A//  ( ka )  4^  4 

Ps^a 


Pa 


cos  (fia2a)  COS  (Qszs)  +  sin  (fia^a)  sin  (fis2s 


=  0 ,m  =  0, 1,  2,  3... 

Change  of  variables  for  analysis  of  the  characteristic  equation  in  MATLAB™  requires  that  rmn  = 
za \Jk‘l  —  (mn,  which  implies  that  £mn  =  k2a  —  ■  This  ensures  the  eigenvalues,  rmn,  will 

be  approximately  evenly  spaced,  similar  to  the  change  of  variables  from  the  two-layer  waveguide. 
The  characteristic  equation  therefore  becomes 


Ps“ 


t  (fc2-U4N 


(8.9) 


x 


PsT 


s  1  mn 


■■  cos  (rmn)  sin  (zsy/k2s  -  Cmn)  -  sin  (rmn)  cos  (z8y/kj  -  ( mn) 


ZaPa\/k2s  ~  C 
X  ^ Jm  ^y/ C ninP^j  4  Jm  {kok)  1 2^  T"  4 

PsTmn  COS  (rmn)  COS  -  Cmn)  +  V' ^  “  Cmn  SHI  (rmn)  sill  -  Cmn) 


X 


zap, 

=  0,  m  =  0, 1,  2,  3... 


8.1.2  Predicting  Resonances 

Consider  equation  (8.9).  Substitution  of  Q.s  results  in 

psu;2  1 


r  (fc2-U4(N 


X 


cos  (rmn)  sin  (zsQs)  -  sin  (rmn)  cos  (zsfis) 


AaPa^s 

x  ^ y/ C rai®j  4  Jm  (ko  j  4^  T  4 


X 


ZaP, 

=  0,  m  =  0, 1,  2, 3. 


COS  (t mn)  cos  (zsf)s)  +  y/ A;2  -  Cmri  sin  (rmn)  sin  (zsQs) 
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For  a  Z-resonance  (r  — >  0),  assuming  pa  <C  ps 


Ps^ 


2  1 


sin  (zs\fk2s  -  k2 )  (J0  (kaa)  h  -  J0  ( ka )  I2) 


T  ( k2-k2a ) 

+J0  ( ka )  /2  (cos  (zs^Jk2s  -  fc2)  )  =  0 


with 


and 


/i  = 


k2  -  k 2 


(kaJ0  (ka)  Ji  (kaa)  -  kJx  (ka)  J0  ( kaa )) 


a 

a  /  9  /  7  N  9  /  7  N  \  a 


h  =  -jr-  (J02  +  Ji  {Ka))  -  t~  (mJi  (kaa)  J0  ( kaa )) . 

L  Ka 


For  an  E-resonance  ((mn  — > ►  0), 

,2  x  Ji  (Zca 


T 


(sin  (zsks)) 


a  \  ,  ,  2  -  ' '  '  a 


^  -  Jo  (fca)  -  J  —  /cs/rJ0  (fca)  -  cos  (zsks)  =  0. 


Recall  that  r  =  za^k2a  —  ( mn ,  and  ('  represents  the  resonant  properties  of  the  membrane.  These 
equations  are  plotted  over  a  broad  range  of  frequencies  for  the  Sabatier  parameters  (fig  (8.1)). 


8.1.3  The  Radial  Continuity  Conditions 

The  first  of  the  radial  continuity  conditions  (8.6)  gives 

OO  OO 

E  E  Hm  {\rK°)  ( C'mn  cos  (m6)  +  Dmn  sin  (m0))  Zn  (z) 


m= 0 n= 1 

OO  OO 


EE-; 


(v&o)  Jm  (v7Ea)  COS  (m  (60  -  0)) 


2=0  n=  1  2 £m  ^ r 

xx(z0)  Zn(z0)  Zn(z) 

OO  OO  /  _  \ 

=  Jm  (  )  (Amn  cos  (ra<9)  +  Bmn  sin  (ra<9))  0mn  (z) 


m= 1 n=l 


where 


z  = 


sin  (Oa  (z-  za)) 

COS  (ELaza)  sin  (Osz)  -  sin  (Qa  (za))  cos  (O. 


0  <  z  <  za 
\sz)  zs  <  z  <  0 


The  second  of  the  radial  continuity  conditions  (8.6)  gives 

OO  OO 

E  E  {\JKia)  ( Cmn  cos  (md)  +  Dmn  sin  (md))  Zn  (z) 


m= 0 n=  1 

oo  oo  i  /  n 

E  E  -tNE11™  (VKra)  j'm  (V*v»)  cos  (rn  (»o  -  «))  X  (so)  Z„  (so)  Z„  (z) 


J— 0  71=  1  2<Sm<J) 

OO  OO  _  /  _  \ 

=  E  E  Vc^nj'm  (  VCBa)  (4  cos  (m0)  +  Bmn  sin  (m0))  (pmn  (z) . 

772=1  n=  1  '  ' 


Eigenfunction  Response 
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K  iq'3  x  Eigenvalues  for  the  Sabatier  Parameters 


Frequency  (Hz) 


Figure  8.1:  Eigenvalues  of  the  Membrane  Problem.  The  eigenvalues  corresponding  to  frequencies 
are  determined  from  the  zero  crossings. 
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Equating  the  coefficients  of  the  trigonometric  functions  gives 


uu  /  _  \ 

Jm  y/ Cmn^  )  0mn  M  Amn  -  H «  (VKa)  Z„  (z)  C„ 


2c„/h 


(v7V"o)  (yTV*)  cos  C"^o)  X  (z0)  Zn  (A))  Zn  (2) 


UVJ  /  _  \ 

Jm  \/ Cmn®  )  0mn  (z)  -  H<B  yjr„a)  Zn  (z)  D„ 


2em$. 


-^i1}  (v7V"o)  Jm  (a/ZV*)  sin  (™#o)  X  (20)  Zn  (20)  Z„  (*0 


\/ C mn ^ m  Cran^J  4mn  (z)  Amn  -  '/TnH[E  (VKa)  zn  (z)  a, 

-„/X  -ffm1  (VKr<>)  Jm  (VKa)  COS  (m*Cl)  X  (*>)  Z„  (z0)  Z„  (z) 


V Cmn^m  ^ Cran^J  4 mn  (^0  -^mn  VKHm'  {VKa)  ^ n  (^)  Drnn 

=  -Pert  H™  (VKro)  J'm  {VJNa)  sin  (m0o)  X  (-0)  Zn  (*o)  Z„  (*)  . 

If  60  —  0  it  suffices  to  consider  the  cosine  equations.  Multiplying  by  ;\  Z„  and  integrating  on  z  from 
z.5  to  za  and  using  the  orthogonality  of  the  Zn  gives 

E4n(vM  (f  n  H^nl  ^  Tn®)  77 'mn 


=  ~$r 


2£m$. 


-Hm  (v7hiro)  Jm  {ViNa)  x  (zo)  Zn  (z0) 


where,  according  to  MAPLE™, 


E  y/ClJm  y\/Cla)  I'rnnlAml  *&n\/ En^ni  (\/ Ena)  ^’n 
H™  (V^r°)  (v7^0)  X  (*o)  Zn  (z0) 


$„  =  f  x  (z)  Zn  (z)2  dz  —  f  paZSt1l  (z)2  dz  +  f  psZan  (z)2  dz 


-f1  sin(Cn)2(cos(«n)  sin (an)  +  an)  +  cos(an)2(-  cos(fn)  sin(f„)  +  £ 
2  a  2^n 
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‘-mnl 


where 


and 


U  Za 

fx  {z)  4* ml  {z)  Zn  (z)  dz  =  JpJm^s  ( Z )  Zn,s  (z)  dz  +  f ps<t>ml,a  (-)  Zn,a  {z)  dz 

Zs  Zs  0 

sin  (£J  f  -zspsTmn  cos  (. tmn )  cos  (rmn)  +  zapaoin  sin  (an)  sin  (r mn) 
Za  (- i2mn  +  a2n)  V  -ImnZaPa  sin  (w)  sin  (rmn)  +  ZspsTmn  COS  (an)  COS  (t mn) 
zaps  COS  (an)  {—Tmn  cos  (rmn)  sin  (£n)  +  fn  sin  (rmn)  cos  (£J) 

e  -t2 

Sn  mn 


Za^fk^a  /A 

Oln  Zs  \J k2  pr 

Tmn  Za\Jk “  C71 

L’mn  \/  kg  6// 


dn  =  K-['^ 

Zn 


Kr 


Eliminating  Cmri  gives 

°o  r  / 

E  VKHIP'(VK°)  ’km  (  \lClrna 


1=1 


ini  ,Jm 


Im ®>)  ( V Ten)  Imnl-A-{$AIT) 


-AjS-xl*,)  Zn(zo)  H™  (^JTnra) 

ZSm, 


H, 


(1)' 


(v7 Erfl)  km  (v^n°) 
~J'm  (v^0)  ^  (v7/7^) 


Equation  (8.10)  is  an  infinite  dimensional  linear  system  for  Ami .  In  computing  /lm/ ,  the  infinite 
limit  will  be  replaced  by  a  sufficiently  high  value,  N,  to  fonn  a  finite  dimensional  system.  This  is 
one  of  the  equations  that  will  be  programmed  into  MATLAB™  to  solve  the  Membrane  Problem. 
Note  that  C'v  {z)  =  —  Cv+i  (z)  +  ^Cv  {z),  if  C  is  a  Bessel  or  Hankel  Function.2  For  simplicity, 
rewrite  equation  (8.10),  defining  the  first  and  second  lines. 


2[2],pg.  361. 
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Combining  like  tenns  yields 


The  bracketed  portion  of  line  (8.11)  can  be  rewritten 


Combining  like  tenns  yields 


X  -y/KHml  1  (vKa)  - 

y  V  Tna 

Jr\fClrn,Jm+ 1  ^  \/ Clm ^-rn  (V TCa) 

\/ Clm^rn-Y- 1  "s/Clrn^)  ^rn  (\//VT)  Tn  '4/  (  \i Cirri' 


-Hg> 


a)  #1+1  (V/x^a) 


(8.14) 


Lines  (8.13)  and  (8.14)  were  entered  into  the  formula  for  /lmn .  Therefore, 

X/  \Z~Clm.Jm+ 1  ^  \/  Clm0^  #1^  (\/  Prfl)  ~  V  Pn  Jm  ^\/  C  #m+l  (V  Cn°)  Irnnl-A-ml 

=  -+E  U)  Z»  (So)  '!>  (VS’-o)  f  (v/C“(  . 

2£m  Jm  (^/ Ena)  #m+l  (V Ena) 

The  C'mn  are  then  given  by 


G,nii 


oo 

£  -4,  h/Q»)  .g(i)  Jm  (+^a)  x  (.„)  Z„  (+ 

4?n#m  (^/ /Xno)  2  Em^nHm.  (^/  TnP) 


Velocity  in  the  atmosphere  is 


w  =  - chp 

^Pa 

|  OO  OO  /  _  \ 

=  T -  E  E  (  yEmr  )  4n  COS  (m0)  (f)'mn  (-)  ,  T  <  a 

Pa  m=0  n=l  '  ' 


(8.15) 


CHAPTER  8.  MEMBRANE  IMBEDDED  IN  A  RIGID  SUBSTRATE 


63 


where 


0mn  {z)  —  sin  (Oa  (z-za)). 


<\>mn  (Z)  =  COS  (Oa  (z  -  Za))  . 


Chapter  9 


Membrane  Analysis  and  Comparison  with 
Sabatier1 


9.1  Membrane  Analysis 

In  this  section,  the  blue  line  in  each  graph  represents  the  Sabatier  parameters.  All  other  colors 
are  variations  on  the  indicated  parameter.  Note  that  the  Sabatier  parameters  are  zs  =  —0.175  m, 

ca  =  340  m/s,  pa  =  1.2  H,  ccomPjS  =  160  m/s,  cshear,s  =  75  m/s,  and  ps  =  1400  H. 


9.1.1  Termination  of  the  Infinite  Sum 

In  order  to  obtain  reasonable  solutions,  it  is  necessary  to  extend  the  solution  in  equation  (8.15)  to 
at  least  the  first  soil  eigenvalue.  This  eigenvalue  results  from  equation  (4.18),  particularly  from 

the  cos  zs\Jk2  —  p  in  the  first  line.  Since  cos  (x)  =  0  when  x  =  (2"2ll7r,  zsy/k 2  —  p  =  t2"21^. 

Therefore,  p  —  k2s  —  ( K  —  zaJkl~k2+  .  Since  the  first  soil  resonance 


2  zs 


2  zs 


is  of  interest,  n  —  1,  and  =  za  \  /  k2  —  k2  +  (  yr  )  •  With  the  appropriate  substitutions  made, 


2  { 


+  (  ^  )  .  To  ensure  that  the  first  resonance  is  completely 


K  =  5oo,/(^y 

accounted  for,  k  +  1000  will  be  used  to  terminate  the  infinite  sum.  Figure  (9.1)  represents  early 
and  late  termination  of  the  sum,  around  this  predicted  value,  for  the  first  mechanical  resonance  for 
the  Sabatier  parameters. 
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Soil  Particle  Velocity  (m/s) 


CHAPTER  9.  MEMBRANE  ANALYSIS  AND  COMPARISON  WITH  SABATIER2 


Soil  SurfgeftjProfile  using  the  Sabatier  Parameters  with  Early  Termination  of  the  Infinite  Sum 
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Figure  9.1:  This  figure  represents  the  effects  of  early  termination  of  the  infinite  sum.  Note  the 
exaggeration  and  unreasonable  difference  in  magnitude  between  the  resonance.  The  vertical  scale 
is  1  x  10-3. 
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S#i^|§>?jrface  Vibration  Profile  using  Sabatier's  Parameters  with  Depth  Variance 
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Figure  9.2:  This  graphic  shows  how  increasing  the  depth  results  in  much  weaker  mechanical  and 
acoustic  resonances.  Note  how  the  resonances  change  with  depth.  The  vertical  scale  is  1  x  10~6. 
The  units  of  zs  are  meters. 


67 


CHAPTER  9.  MEMBRANE  ANALYSIS  AND  COMPARISON  WITH  SABATIER4 


Figure  9.3:  Results  of  varying  the  compressional  wave  speed  within  the  Sabatier  parameters.  The 
units  of  cPjS  are  m/s. 


9.1.2  Variance  of  Parameters 

Depth  of  Mine  Detection  Limitations 

As  seen  in  Figure  (9.2),  as  depth  increases,  the  mechanical  resonance  moves  to  the  right  and 
decreases  in  strength.  If  the  large  mechanical  resonance  disappears,  the  acoustical  resonance  is 
also  very  hard  to  detect.  The  resonance  of  the  mine  itself  is  much  stronger  than  the  acoustic 
resonance.  Resonances  at  larger  depths  are  soil  resonances,  and  can  be  predicted  from  the  two- 
layer  waveguide.  According  to  Sabatier,  landmines  are  usually  buried  less  than  10  centimeters  (.  1 
m)  from  the  soil  surface.  As  indicated  by  Figure  (9.2),  the  soil  particle  velocity  plots  are  more 
precise  at  shallow  depths. 

Soil  Variations 

Compressional  sound  speed  (fig.  (9.3))  affects  the  curvature,  location,  and  magnitude  of  the  reso¬ 
nance.  Larger  values  for  compressional  sound  speed  move  the  resonances  to  the  left,  and  decrease 
the  curvature  of  the  leading  part  of  the  resonance.  There  does  not  appear  to  be  any  correlation 

'[27] 
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Soil  Si#fflg"<§  Velocity  Plot  for  the  Sabatier  Parameters  with  Shear  Wave  Speed  Variation 


Frequency  (Hz) 


Figure  9.4:  Results  of  varying  the  shear  wave  speed  within  the  Sabatier  parameters.  The  units  of 
cs  are  m/s.  The  vertical  scale  is  1  x  10~5. 


between  large  and  small  values  with  magnitude;  however,  magnitude  does  show  some  dependence 
on  compressional  wave  speed.  Note  that  Sabatier  suggests  a  value  of  160  — ,  while  experiments 
from  Dr.  D.  Keith  Wilson  and  Dr.  Flarley  Cudney  of  the  US  Army  Corps  of  Engineers  suggest 
values  may  be  as  high  as  647  ™  in  shallow  desert  soil.  With  such  a  large  variance  in  soil  parame¬ 
ters,  ascertaining  the  effect  of  sound  speed  on  the  profiles  is  a  very  important  step  in  predicting  the 
landmine  resonances. 

According  to  Figure  (9.4),  larger  values  of  shear  sound  speed  appear  to  move  the  resonance  to 
the  right,  and  decrease  the  intensity  of  the  soil  velocity  response. 

Similarly,  the  density  of  the  soil  appears  to  have  significant  effect  (fig.  (9.5)).  This  parameter 
seems  to  adjust  resonance  width,  amplitude,  and  a  slight  shift  in  the  resonant  frequency.  The  only 
apparent  correlation  is  between  the  magnitude  of  pp  and  the  resonance  shift:  small  pp  appear  to 
the  right  of  larger  values  for  pp.  Sabatier  suggests  a  value  of  1 400 while  the  experimental  data 
from  Wilson  and  Cudney  suggest  1600  Note  that  ps  =  pp  +  (A  and  that  the  pp  was  altered 
by  changing  ps;  the  value  given  the  key  for  the  graph. 


69 


CHAPTER  9.  MEMBRANE  ANALYSIS  AND  COMPARISON  WITH  SABATIER6 


^flif^urface  Velocity  Profile  for  the  Sabatier  Parameters  with  pp  Variance 


Frequency  (Hz) 

Figure  9.5:  This  graphic  demonstrates  the  sensitivity  of  the  soil  particle  velocity  to  changes  in 
pp.The  vertical  scale  is  1  x  10“6. 
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SojJ  iffyfrace  Velocity  Profile  for  the  Sabatier  Parameters  with  Varying  Mine  Size 


Frequency  (Hz) 

Figure  9.6:  This  profile  shows  the  effect  of  the  size  of  the  mine  with  regard  to  the  resonance.  Note 
even  small  changes  in  the  radius  a  seem  to  have  large  effects.  The  vertical  scale  is  1  x  10~5. 


9.1.3  Membrane  Variation 

The  membrane  density  and  tension  are  difficult  quantities  to  assign  since  they  are  particular  to 
each  landmine.  Since  this  model  is  based  on  the  VS  1.6  Anti-Tank  mine,  a  brief  discussion  of 
the  selection  process  for  pm  follows  in  the  discussion  on  the  comparison  with  Sabatier’s  model. 
However,  this  discussion  suggests  that  either  the  mechanical  resonance  or  exact  parameters  with 
regard  to  both  tension  and  density  must  be  known  for  accurate  detection. 

The  size  of  the  landmine  also  affects  detectability.  The  VS  1.6  is  approximately  the  size  of  a 
dinner  plate,  with  a  radius  of  7  centimeters  (.07  m).  However,  the  smaller  the  mine,  the  smaller 
the  amplitude  of  the  resonances,  making  the  mine  harder  to  detect.  This  is  shown  in  Figure  (9.6). 
Note  that  small  changes  (not  even  order  of  magnitude  changes)  have  a  large  effect  on  the  resonant 
frequencies. 


71 


CHAPTER  9.  MEMBRANE  ANALYSIS  AND  COMPARISON  WITH  SABATIERU 

9.2  Comparison  with  Sabatier’s  Findings 

9.2.1  The  Sabatier  -  Waxier  -  Velea  Paper7 

In  An  Effective  Fluid  Model  for  Landmine  Detection  using  Acoustic  to  Seismic  Coupling,  Sabatier 
et.  al  create  a  "right  circular  waveguide  with  rigid  walls,  [which]  contains  air  in  the  upper  half, 
soil  in  the  lower  half,  and  a  buried  mine  placed  concentrically  on  the  waveguide’s  axis."8  Using 
a  mass-spring  problem  analogy,  Sabatier  then  uses  a  "brute  force  numerical  scheme"9  to  solve  the 
pressure  equation,  using  the  same  boundary  conditions  as  this  project.  The  paper  from  Sabatier  et 
al.  does  not  consider  the  scattering  effects  of  the  mine  or  soil  for  an  incident  acoustic  or  seismic 
wave;  whereas  this  project  has  included  the  scattering  effects. 

The  Sabatier  parameters  refer  to  the  properties  used  by  Sabatier  in  his  research.  The  first 
parameter  described  is  the  depth  of  the  buried  mine,  referred  to  as  zs  in  this  project,  at  3  inches,  or 
0.075  meters.  The  mine’s  radius  is  next  mentioned,  referred  to  as  a  in  this  project,  and  assigned 
a  value  of  0.07  meters.  The  frequency  sweep  that  Sabatier  perfonns  ranges  from  /  =  80  —  300 
Hz.  Finally,  the  air  and  soil  parameters  that  Sabatier  uses  are  listed  as  ca  =  340  m/s,  pa  =  1.2  ^4, 
cs  =  160  m/s,  and  ps  =  1400 

9.2.2  Selection  of  pm,  the  Density  of  the  Membrane 

The  selection  of  prn,  along  with  other  variables,  detennines  the  resonances  for  the  soil-membrane 
system.  Sabatier’s  model  of  the  system  as  a  mass-spring  damper  lists  parameters  specific  to  the 
mass  and  spring  that  cannot  be  easily  translated  into  values  for  this  project.  Expecting  that  the 
mechanical  resonance  should  occur  at  approximately  100  Hz,  reasonable  approximations  are  made 
from  the  mathematical  equations  for  pm.  From  equation  (8.9),  a  mechanical  resonance  occurs 

when  J0  ( ka )  =  0,  where  feel.  Recall  that  k  =  js  actually  complex,  where  u  is 

the  angular  frequency  (u  =  2nf),  [3rn  accounts  for  the  damping,  and  T  is  the  tension.  The  first 
zero  of  the  Bessel  function  occurs  when  ka  =  2.405. 11  If  a  =  .07  m,  as  in  Sabatier’s  experiment, 

then  k  =  34.3571  «  \j!>"f  .  Sabatier’s  parameters  are  assumed  for  the  other  variables  to  be 
u  =  2vr  (100)  Hz,  (3m  =  1000  N-s/m,  and  T  =  2,  000,  000  kg/s2,  then  pm  =  5980.04  «  6000 
kg/m3. 

9.2.3  Comparison  with  Sabatier’s  Work 

A  mechanical  resonance  is  defined  as  "the  resonant  frequencies  of  any  mechanical  system  for 
which  the  input  mechanical  reactance  goes  to  zero."12  The  mechanical  system  refers  to  the  land¬ 
mine  itself,  explaining  the  narrow  band  response  observed  in  this  paper’s  frequency  response  graph. 

7[27] 

8[27],  pg.  1993. 

9[27],  pg.  1996. 

11  [5],  pg.  x. 

12[18],  pg.  48. 


Velocity  magnitude  (  m/s ) 
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Frequency  (  Hz ) 


Figure  9.7:  Sabatier’s  measured  values  in  a  soil  surface  velocity  plot.  Note  the  increments  are 
approximately  10  Hz. 


Soil  Particle  Velocity  (m/s) 
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Soil  Surface  Velocity  Profile  for  the  Sabatier  Parameters 


Figure  9.8:  Predictions  based  on  the  theoretical  solution  presented  in  this  paper  for  Sabatier’s 
parameters. 
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Man-made  objects,  like  a  landmine  case,  tend  to  have  a  relatively  narrow  band  response  due  to  the 
highly  organized  molecular  structure  of  the  solid  material  when  compared  with  the  granular  and 
significantly  more  random  structure  of  the  soil  particles.  This  higher  quality  factor  "Q"  response 
allows  acoustic  to  seismic  landmine  detection  to  remain  fairly  accurate  (in  discriminating  detec¬ 
tions  from  false  alanns)  in  even  rocky  soil,  where  conventional  detection  techniques  fail.  Note  the 
common  resonance  in  the  Sabatier  plot  (fig.  (9.7))  and  this  project’s  plot  (fig.  (9.8))  at  approxi¬ 
mately  100  Hz. 

The  anti-resonance  is  the  frequencies  which  result  in  a  local  minimum  for  the  frequency  re¬ 
sponse  plot.  This  can  be  compared  to  the  destructive  interference  between  two  propagating  waves. 
Sabatier  claims  the  anti-resonance  occurs  from  1 80-2 1 0  Hz,  while  the  data  from  this  paper  suggests 
anti-resonances  occur  precisely  before  the  mechanical  resonances. 

The  prediction  also  suggests  a  second  mechanical  resonances  for  the  membrane  at  approxi¬ 
mately  230  Hz. 

There  appears  to  be  a  peak  in  Sabatier’s  work  near  /  =  240  Hz  that  is  unlabeled.  This  is  a  soil 
resonance,  predicted  by  the  two-layer  waveguide  problem  in  this  paper  (see  fig.  (9.9)).  The  same 
soil  resonance  is  observed  in  the  plot  from  this  project’s  derivation. 

The  first  acoustical  resonance  is  tenn  defined  by  Sabatier  to  be  the  frequency  at  which  "the  layer 
of  soil  between  the  top  of  the  mine  and  the  surface  of  the  soil"16  vibrates  at  maximum  amplitude. 
This  resonance  occurs  very  clearly  in  Sabatier’s  graph  at  approximately  /  =  280  Hz.  However, 
according  to  Figure  (9.9),  the  second  soil  resonance  occurs  at  this  same  frequency. 


16[27],  pg.  1994. 


Soil  Particle  Velocity  (m/s) 
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x  iq'8  False  Alarm  Plot  for  the  Sabatier  Parameters 


Figure  9.9:  This  plot  shows  the  predicted  soil  resonances  for  the  false  alarm  case,  generated  from 
the  two-layer  waveguide  using  the  Sabatier  parameters. 


Chapter  10 

COMSOL™  Advances 


COMSOL™  is  a  sophisticated  multiphysics  program  which  employs  the  finite  element  method  to 
solve  very  complicated  problems  of  multiple  equations  very  quickly.  The  finite  element  method 
discretize  the  domain  into  smaller  areas,  and  solves  the  equations  over  these  small  areas.  For 
clarification,  an  example  is  included  in  Figure  (10.1).  COMSOL™  is  well-known  for  its  selec¬ 
tive  meshing  technique,  which  pennits  the  user  to  use  rough  meshing  in  less  important  parts  of 
the  diagram,  and  fine  meshing  in  important  parts  of  the  diagram.  The  diagram  used  to  calcu¬ 
late  these  results,  along  with  its  selective  meshing,  is  shown  in  Figure  (10.1).  For  the  two-layer 
waveguide,  the  two-dimensional  axially  symmetric,  time  hannonic  pressure  acoustics  module  is 
used.  The  Bessel  panel  and  loudspeaker  tutorials  suggested  the  spherical  geometry  developed  in 
Figure  (10.2)  to  represent  the  two-layer  waveguide.  The  Bessel  panel  tutorial  was  chosen  because 
of  its  implementation  of  the  delta  function.From  the  data  obtained  in  the  solution,  shown  as  the 
subplot  within  Figure  (10.2),  multiple  post-processing  options  are  available.  Of  particular  interest 
for  this  problem  are  the  radial  distance  versus  pressure  plots,  an  example  of  which  is  shown  in 
Figure  (10.3). 

Unfortunately,  COMSOL™  has  not  produced  results  that  agree  with  the  two-layer  waveguide 
results.  The  main  problem  appears  to  be  the  implementation  of  the  delta  function  in  COMSOL™. 


Figure  10.1:  The  illustration  above  demonstrates  how  the  Finite  Element  Method  would  calculate 
the  perimeter  of  a  circle.  The  middle  circle  would  be  considered  a  rough  mesh,  while  the  circle  on 
the  right  would  give  a  more  accurate  solution. 
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Figure  10.2:  This  diagram  shows  the  geometry  of  the  COMSOL  model  drafted  for  the  two-layer 
waveguide.  A  quarter  circle  with  a  radius  of  500  representing  the  atmospheric  layer  was  placed 
atop  a  rectangle  of  depth  1  meter.  Appropriate  boundary  conditions  were  selected,  and  the  solution 
in  the  location  around  the  point  source  is  shown  as  an  insert. 


The  source  in  the  problem  plotted  in  Figure  (10.3)  was  at  r0  =  0,  so  the  highest  pressure  should 
logically  appear  at  r  ~  0.  The  COMSOL™  solution  is  still  being  worked  on,  and  hopefully  a 
result  matching  the  MATLAB™  computed  two-layer  waveguide  will  be  obtained. 
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Acoustic  Pressure 


Figure  10.3:  An  example  of  a  radial  distance  versus  pressure  plot  for  the  two-layer  waveguide. 
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Appendix  A 

Orthogonality  Proofs 


A.l  Orthogonality  of  the  Eigenfunctions 

The  generic  Sturm-Liouville  Differential  Equation  in  the  z  coordinate  is  represented  by 

(p(z)  +  [**(*)  +  q(z)\z  =  °-  (A.l) 

In  the  case  of  the  two  layer  waveguide  problem, 

z;(z)  +  (-»  +  k*)Z'(z)  =  o 
Z'J(z)  +  (-^  +  k2,)Z,(z)=  0 

p(z)  —  1,  A  =  —  p,  and  q  (z)  —  |  ^  ^  ^ z^<{)  ’ t0  °^ta^n  ^  (z 

function.  Let  Zm  and  Zn  be  eigenfunctions  of  the  depth  equation  (Z  (z)).  The  following  proof 
implies  that  eigenfunctions  Zm  and  Zn  of  the  Sturm-Liouville  equation  are  orthogonal  on  (a,  b ) 
with  respect  to  the  weight  function  x  (z)  if 

b  b  b 

J Zm  (z)  ( p  (z)  Z'n  (z))'  dz  +  fq  (z)  Zm  (z)  Zn  (z)  dz  +  XnJ x  {z)  Zm  (z)  Zn  (z)  dz  =  0.  (A.2) 

a  a  a 

In  the  case  of  the  two  layer  waveguide  problem,  the  condition  above  refers  to  the  two  intervals  of 
(zs,  0)  and  (0,  za),  written  as 

oo  o 

f  ZS)m  (z)  Z"s  n  (z)  dz  +  f  k2ZSt,n  (z)  ZS)H  (z)  dz  -  pnJ x  {z)  ZS:Tn  (z)  Zs,n  (z)  dz  =  0 

£s  Zs 

Za  Za  Za 

f  Zam  (z)  Z"a  n  (z)  dz  +  J  k2Zarn  (z)  Za,n  (z)  dz  -  pnf  x  (z)  Za,m  (z)  Za,n  (z)  dz  =  0. 

0  ’  0  0 


(  Za  0  <  z  <  za 
\  Zs  zs  <  z  <  0 
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Using  integration  by  parts  on  the  first  integral  of  both  conditions  gives 

JZm  (z)  (P  (z)  Z'n  (z))'  dz  =  Zm  (z)  p  (z)  Z'n  (z)  \ba  -  fp  (z)  Z'm  (z)  Z'n  (z)  dz 

a  a 

or,  in  the  specific  cases  presented  in  this  paper, 

fZs,m  ( z )  Zs  n  (z)  dz  =  ZSyTn  (z)  Zs  n  (z)  \Zs  —  f  Zs  m  (z)  Zs  n  (z)  dz 


I Za,m  M  Z"„  W  dz  =  Za.m  (z)  Z’a  n  (z)  -  f  Z'm  (z)  z;„  (z)  dz. 

0  0 

The  governing  boundary  conditions  are  that  Z  (za)  =  0,  Za  (0)  =  Zs  (0),  psZ'a  (0)  =  paZ's  (0)  and 
Z'  (zs)  =  0.  Therefore, 

Zs,m  {z)  Zsn  (z)  \Zg  =  ZSyTn  (0)  Zs  n  (0)  —  ZS)m  (zs)  Zs  n  ( zs )  =  Zsm  (0)  Zs  n  (0) 

Za,m  {z)  Z'a  n  (z)  1 5“  =  ZayTn  (za)  Z'a  n  (za)  —  ZayTn  (0)  Z'an  (0)  =  —  za,m  (0)  Z'a  n  (0) 

=  (0)  z;,„  (0) . 

Ps 


Z s,m  (o)  Z'sn  (0 )-  fZ'sm  (z)  Z'sn  (z)  dz+fk2ZSym  (z)  ZS)n  (z)  dz-pn  jx  (-)  Zs,m  {z)  ZSyU  (z)  dz  =  0 


0  —  Za,m  (0)  Za  n  (0)  JZa>m  (z)  Za  n  (z)  dz 

o 


T  f  kaZarn  (^)  ^ a,n  (^)  dz  PnJ X  (^)  ^ a,m  (^)  ^ a,n  (^)  dz. 

0  0 


Switching  m  and  n  gives 


0  —  ZSyH  (0)  Zs  m  (0)  J Zs  n  (z)  Zs  m  (z)  dz 

Zs 

0  0 
T  J^s^s,n  iz)  ^ s,m  (^)  dz  PmJ X  (^)  ^ s,n  (^)  ^s,m  (^)  dz 


=  -zayn  (0)  Z'a  m  (0)  -  f  Z'a  n  (z)  Z'a  m  (z)  dz 


'■'a  '■'a 

“I”  f  (^0  ^ a,m  (^)  dz  l^mf  X  (^)  ^a,n  (^)  ^ a,m  (^)  dz. 
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Subtracting  the  two  equations,  respectively,  results  in 


Zs,m  (o)  Z's  n  (0)  -  Zs,n  (0)  Z's  m  (0)  -  (nn  -  fim)  fx  (z)  Zs>n  (z)  Zs 


(z)  dz  =  0 


-Za,m  (0)  z'an  (0)  +  zain  (0)  Z'am  (0)  -  {pn  -  nm)  fx  (z)  Za,n  (z)  Za?m  (z)  dz  =  0. 

0 

Note  that  the  second  equation  above  is  actually 

[Zs,m  (0)  Z'  (0)  -  (0)  Z'  (0)]  -  (pn  -  pj  fx  (z)  Za,n  (z)  Zam  (z)  dz  =  0 

Ps  0 

or 

-Zs,m  (o)  z'sn  (0)  +  ZS)U  (0)  Z'sm  (0)  -  —  pm)  fx  (z)  Za,n  ( z )  Za,m  (z)  dz  =  0 

Pa  0 

so  that  the  first  equation  added  to  the  equivalent  second  equation  is 

0  P  Za 

—  {Pn  Pm)  f  X  Zs^n  (z)  Zsrn  (z)  dz  —  (/rn  Pm)  f  X  {z)  Za>n  (z)  Za^n  (z)  dz  =  0 

zs  Pa  0 

0  Za 

-Pa  ( Pn  -  Pm)  fx  (-)  ZS)U  {z)  Zs,m  (z)  dz  ~  ps  (pn  -  flm)  fx  {z)  Za,n  (z)  Za>m  (z)  dz  =  0 


( Pn  Pm) 


Paf  X  (z)  ZS)U  (z)  Zsm  (z)  dz  +  pjx  (z)  Za)U  (z)  Zam  {z)  dz 

Zs  0 


=  0. 


Assuming  \n  f  Xm  when  n  f  m. 


fx  ( z )  Zn  (z)  Zm  (z)  =  0 

Zs 


where 


Zn  (z)  Zm  (z) 


Za,n  (z)  Za  m  (z)  dz  when  0  <  z  < 

Zs  n  (z)  Zs^m  (z)  dz  when  zs  <  z  <  0 


implies  Z  is  orthogonal  with  respect  to  the  weight  function 


X(z) 


ps  when  0  <  z  <  za 
pa  when  zs  <  z  <  0 


Appendix  B 
Additional  Proofs 

B.l  Explanation  of  the  Exclusion  of  Gravity  from  the  Conser¬ 
vation  of  Momentum  Equation 


The  exclusion  of  the  gravity  tenn  from  the  derivation  of  momentum  originates  from  the  following 

Although  gravity  is  always  present,  it  has  negligible  influence  on  acoustic  distur¬ 
bances  of  all  but  extremely  low  frequencies,  e.g.,  those  of  order  or  less  than  where 
g  is  the  acceleration  due  to  gravity  and  c  is  the  speed  of  sound;  so,  for  simplicity,  the 
body  force  tenn  is  neglected  at  the  outset.  Acoustic-gravity  waves  [internal  waves] 
(infrasonic  waves  with  frequencies  so  low  as  to  be  strongly  affected  by  gravity)  have 
been  a  major  topic  of  research  during  the  past  two  decades  but  fall  outside  the  scope 
of  an  introductory  discussion.1 


Since  cs  =  160  m/  s,  and  g  =  9.8  m/  s2,  /  =  9  g°^/s2  =  16.  327  Hz  100  Hz,  which  is  the 
lowest  frequency  of  interest. 


B.2  Proof  of  the  Denominator  in  equation  (4.25)  Using  Hand¬ 
book  of  Mathematical  Functions 2 

The  denominator  first  reads 

(v'/V'o)  Jk  (vTfl^o)  +  VFiHjp  (V/fiA))  J'k  (V/flr o) 

By  Abramowitz  and  Stegun,3  the  derivative  of  a  Bessel  or  Hankel  function  can  be  simplified  into 


‘[4],  [13],  [15],  [10] 

2[2],pg.  362. 

3[2],  pg.  362,  eqn.  (9.1.27). 
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-vs ( ~ OENNi PEA  +  (v^r°) )  A (v^ro) 

V  v^r  0  / 

-  (v^t-q)  Jfc+1  (y^Ao)  +  kJk  (yfiHro)  \ 
VWir<  o  / 


+v/^i/l1)  (V^o) 


Simplifying  the  denominator  of  the  equation  yields 

(v^ro)  fern  {V^ir o)  (V^ro)  -  ^  (v/^ro)  Jfc+i  (Vfh'ro) 


To 


k 


Hk]  (v/^ro)  {\/THro)  -  Hl1}  (v^qro)  Jfc  (y^o) 


(i) 


+  0 


r0 

Recognizing  that  the  second  part  of  the  equation  is  zero  gives 

(vTh)  Hkl i  (vTh^o)  Jk  (y^o)  - 4  5  (y7¥~o)  Jfc+i  (y7qr0) 

Multiplying  to  use  Abramowitz  and  Stegun,4 

|  #1+1  {\/Jhro)  H[}]  (y^qro)  +  f/fi  (y^A))  i?*0  (y^A)) 
"#*  (v^ro)  #$!  (y^A>)  -  ^  (y^o)  ifg,  (y^o) 


Applying  Abramowitz  and  Stegun5  gives 


\TEi 


Nil  (V/Vo)  N>  (Vw)  -  A”  (VSn>)  A«  (VSro) 


r(2) 


(i) 


r(2) 


2* 


Vri  f 

2  V  ^yW  ^0 


4*  \ 


4[2],  pg.  362,  eqn.  (9.1.17). 

5[2],  pg.  362,  eqn.  (9.1.17). 


Appendix  C 

Graph  Parameters 


C.l  Common  Values 

Note  that  units  have  been  dropped  for  the  numerical  analysis 
=  500 
ca  =  330 
Pa  =  1-168 
a  =  0 

rtol  =  5  x  10~5 
Pm  =  1000 
T  =  2  x  106 
a  =  0.07 
m  =  0 


C.2  Values  for  Specific  Figures 

Figures  (4.2)  and  (8.1): 

=  -0.075 
cs  =  170 
cp  =  75 

To  =  2 
Zq  =  2 

0o  =  O 
r  =  0.0 
z  =  0 
0  =  0 

Figures  (4.3),  (4.4),  (6.1),  (6.2),  (6.3),  (6.5),  (6.6),  (6.7),  and  (6.4): 
=  -1 
cs  =  160 
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T  0  =  2 
zq  =  2 

0O  =  O 
r  =  0.0 
z  =  0 
0  =  0 

Figures  (9.1),  (9.2),  (9.3),  (9.4),  (9.5),  (9.6),  and  (9.8): 

=  -0.175 
cs  =  170 
cp  =  75 
f o  =  0.2 

^o  =  2 

0o  =  O 
r  =  0.0 
z  =  0 
0  =  0 


Appendix  D 
MATLAB™  Files 

D.l  Common  Functions 

D.l.l  csqrt.m 

function  y=csqrt(z) 
y=sqrt(z); 

iset=find(  imag(z)  <  0) ; 
y(iset)=-y(iset); 

D.1.2  eigs4z.m 

function  Z=eigs4z(kap,rhoa,rhos,za,ca,zs,cs, sigma, rhop,f) 
omega=2*pi*f; 

ks=sqrt((omegaA2+i*  omega*  sigma/rhop)/csA2) ; 
ka=omega/ca; 

%ps=-i*omega*rhop/(sigma+i*omega*rhop); 
llamb  =  sqrt(ksA2-kaA2+kap.A2/zaA2); 

Z=rhos*(cos(kap).*cos(zs*llamb)).*(kap/za)+rhoa*(sin(kap).*sin(zs*llamb)).*llamb; 

D.1.3  mu.m 

function  [mu, kappa]=mu(n,ka,rhoa,rhos,za,ca,zs,cs, sigma, rhop,f) 
ftntol=5e-4; 

860101=56-4; 
count=l; 
for  ii=n 

[kap,fv]  =  secant(@(kapl)  eigs4z(kapl,rhoa,rhos,za,ca,zs,cs, sigma, rhop,f), 

(2*ii-l)*pi/2-pi/8,(2*ii-l)*pi/2+pi/8,sectol); 

if  (fv<  ftntol)&(real(kap)  >  (ii- 1 )  *pi)&(real(kap)  <  ii*pi) 
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kappa(count)=kap; 

else 

kl=(2*ii-l)*pi/2-1.57; 
k2=(2*ii-l)*pi/2+1.57; 
kapr=linspace(k  1  ,k2 , 1 00); 

y=real(eigs4z(kapr,rhoa,rhos,za,ca,zs,cs, sigma, rhop,f)); 

if  y(l)>0 

I=find(y<0); 

else 

l=find(y>0); 

end 

if  isempty(I) 

if  count  ==  1 

kappa(count)  =  0; 

else 

Figure 

plot(kapr,y) 

error(’zero  not  found’) 

end 

else 

[kap,fv]=secant(@(kapl)  eigs4z(kapl,rhoa,rhos,za,ca,zs,cs, sigma, rhop,f), 

kapr(I(  1))-.  1  ,kapr(I(  1))+.  1  ,sectol); 

kappa(count)=kap; 

end 

end 

count=count+l; 

end 

if  abs(kappa(l))<sectol 
kappa=kappa(2 :  end) ; 
end 

mu=kaA2  -(kappa/ za) .  A2 ; 

D.1.4  rho.m 

function  wrho=rho(z,rhos,rhoa) 
if  z>=0 
wrho=rhos; 
else 

wrho=rhoa; 

end 
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D.2  Two-Layer  Waveguide  Specific  Programs 

D.2.1  dZn.m 

function  dZ=dZn(z,ka,muv,za,ks,zs) 
if  z>=0 

dZ=(kaA2-muv).A(l/2).*sin((kaA2-muv).A(l/2)*za)./cos((kaA2-muv).A(l/2)*za) 

.*sin((kaA2-muv).A(l/2)*z)+(kaA2-muv).A(l/2).*cos((kaA2-muv).A(l/2)*z); 

else 

dZ=sin((ksA2unuv).A(l/2)*zs)./cos((ksA2-rnuv).A(l/2)*zs).*sin((kaA2-muv).A(l/2)*za) 

./cos((kaA2-muv).A(l/2)*za).Hs((ksA2-muv).A(l/2).*cos((ksA2-muv).A(l/2)*zs) 

./sin((ksA2-muv).A(l/2)*zs).*sin((ksA2-muv).A(l/2)*z)+(ksA2-muv).A(l/2) 

.*cos((ksA2-muv).A(l/2)*z)); 

end 

D.2.2  eigsplot.m 

za=500;  %rough  estimate  of  the  upper  limit  of  the  problem,  units  of  m 
ca=330;  %  speed  of  sound,  m/s 

zs=- 1 ;  %rough  estimate  of  lower  limit  of  problem,  units  of  m 
cs=160;  %Attenborough  et  al.  Attached  derivation  from  Buchanan  units  m/s 
sigma  =  0;  %3 66000; 
rhop  =  4150; 

rhoa=  1.168;  %Wikipedia,  density  of  dry  air  at  std  ambient  press  and  temp 

%Density  of  Air  Article,  units  of  kg/mA3 

r0=2; 

z0=2; 

theta0=0; 

NBlock  =  250; 

r=i; 

theta=0; 

p=[]; 

rtol  =  5e-5; 
f=50; 

omega=2*pi*f; 

rhos=rhop+i*sigma/omega; 

kappa=linspace(0,20,100); 

Y=eigs4z(kappa,rhoa,rhos,za,ca,zs,cs, sigma, rhop, f); 
plot(kappa,Y,’b’) 
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D.2.3  gui.m 

disp(’ Predicted  Soil  Resonances  for  Landmine  Detection’) 

disp(  *******************************************************************  ) 
disp(’ Effective  Fluid  Model  developed  by  Mattingly  and  Buchanan’) 
disp(’  ’) 

disp(’This  method  for  determining  soil  resonances  uses  linearized  Partial’) 
disp(’ Differential  Equations  to  approximate  the  behavior  at  the  air-soil’) 
disp(’ interface.’) 
disp(’  ’) 

disp(  *******************************************************************  ) 
disp(’  ’) 

disp(’All  input  values  are  in  the  assigned  units’) 

disp(’  ’) 

za=500; 

rsp=input(’For  advanced  selections,  press  any  numeric  key  now’); 
if  isempty(rsp) 

typ2=input(’  Enter  the  ambient  temperature  (C)’); 

if  typ2<=-10 

rhoal=1.341; 

cal=325.4; 

elseif  (typ2>-10)&(typ2<=-5) 

rhoal=1.316; 

cal=328.5; 

elseif  (typ2>-5)&(typ2<=0) 

rhoal=1.293; 

cal=331.5; 

elseif  (typ2>0)&(typ2<=5) 
rhoal=1.269; 

031=334.5; 

elseif  (typ2>5)&(typ2<=10) 

rhoal=1.247; 

cal=337.5; 

elseif  (typ2>10)&(typ2<=15) 

rhoal=1.225; 

cal=340.5; 

elseif  (typ2>15)&(typ2<=20) 

rhoal=1.204; 

cal=343.4; 

elseif  (typ2>20)&(typ2<=25) 

rhoal=1.184; 

cal=346.3; 

elseif  (typ2>25)&(typ2<35) 
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cal=349.2; 

rhoal=1.164; 

else 

rhoal=l.l; 

cal=331.4+.6*typ2; 

end 

typ3=input(’Enter  the  humidity  in  decimal  form’); 

rhoa=rhoal*(l+typ3)/(l+1.609*typ3); 

ca=cal+typ3*.6*cal; 

typl=input(’ Select  soil  type:  l=sand,  2=silt,  3=clay,  4=organic  matter’); 
typl  l=input(’Select  soil  moisture  (scale  of  1-5,  l=dry,  5=saturated)’); 
typl2=input(’Select  soil  compaction  (scale  of  1-4,  l=loose,  4=dense)’); 
if  typ  l  ==  1 
if  typ  1 1  ==  1 
cs=260; 

elseif  typ  1 1  ==  2 
if  typ  12  <3 
cs  =  138; 
else 

cs  =  103; 
end 

elseif  typl  1==3 
if  typ  12  <  3 
cs=122; 
else 

cs=253; 

end 

else 

cs=189; 

end 

elseif  typ  1  ==  2 
if  typl  1  ==  1 
if  typ  12  <3 
cs=140; 
else 

cs=139; 

end 

elseif  typ  1 1  ==  2 
if  typ  12  <3 
cs=  117; 
else 

cs  =  122; 
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end 

elseif  typ  1 1  ==3 
if  typ  12  <  3 
cs=121; 
else 

cs=176; 

end 

else 

cs=207; 

end 

elseif  typ  1  ==  3 
if  typ  1 1  ==  1 
if  typ  12  <3 
cs=177; 
else 

cs=159; 

end 

elseif  typ  1 1  ==  2 
if  typ  12  <3 
cs  =  118; 
else 

cs  =  136; 
end 

elseif  typ  1 1  ==3 
if  typ  12  <  3 
cs=107; 
else 

cs=245; 

end 

else 

cs=l  18; 
end 

elseif  typ  1  ==  4 
if  typ  1 1  ==  1 
if  typ  12  <3 
cs=153; 
else 

cs=121; 

end 

elseif  typ  1 1  ==  2 
if  typ  12  <3 
cs  =  89; 
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else 

cs  =  147; 
end 

elseif  typ  1 1  ==3 

if  typ  12  <  3 

cs=87; 

else 

cs=86; 

end 

else 

cs=102; 

end 

else 

disp(’ Error  in  submission  of  soil  type’) 
end 

cp=270; 

sigma=0; 

end 

rtol  =  5e-4; 

NBlock  =  250; 
if  ~isempty(rsp) 

ca=input(’ Enter  the  speed  of  sound  in  the  atmosphere  (m/s)’); 

if  isempty(ca) 

ca=343; 

end 

cp=input( ’Enter  the  speed  of  shear  sound  waves  in  the  soil  (m/s)’); 
if  isempty(cs) 
cs=2 10.31; 
end 

cs=input(’ Enter  the  speed  of  compressional  sound  waves  in  the  soil  (m/s)’); 

if  isempty(cp) 

cp=647; 

end 

sigma=input( ’Enter  the  flow  resistivity  in  the  soil  (Pa*s/mA2)’); 

if  isempty(sigma) 

sigma=0; 

end 

rhoa=input(’ Enter  the  density  of  the  atmosphere  (kg/mA3)’); 
if  isempty(rhoa) 
rhoa=  1.205; 
end 

rtol=input(’ Enter  the  tolerance’); 
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if  isempty(rtol) 

rtol=5e-5; 

end 

end 

NBlock=250; 

zs=input(’  Enter  the  depth  of  the  soil  (m)’); 

if  isempty  (zs) 

zs=-l; 

end 

AA=input(’ Enter  the  position  of  the  source 

[radial  component,  azmuthal  component,  height  component]  (m)’); 

rO=AA(l); 

thetaO=AA(2); 

zO=AA(3); 

r=0; 

theta=0; 

z=0; 

CC=input(’ Enter  the  frequency  range  of  interest 
[lowest  frequency, highest  frequency]  (Hz)’); 

p=ll; 

rhop=8.1827e7/(csA2*.27); 
if  isempty(CC) 

CC=input(’ Enter  the  height  range  of  interest  [lowest  height, highest  height]  (m)’); 
f=input(’ Enter  the  frequency  of  interest  (Hz)’); 
cnt  =  1 ; 

zrange=linspace(CC(l),CC(2),10); 
disp(’  Depth(m)  Pressure  (Pa)’) 
for  z=zrange; 

pp(cnt)=press(p,rtol,NBlock,f, sigma, rhop,cs,ca,rhoa,za,zs,z,zO,r,rO, theta, thetaO); 

disp([z,abs(pp(cnt))]) 

cnt  =  cnt+1; 

end 

plot(zrange,abs(pp),’-m’) 
hold  on 

xlabel(’ Height/Depth  Range  (m)’) 
ylabel(’ Pressure  (Pa)’) 

title(’Soil  Resonance  Predictions  -  Depth  Variation’) 
end 

if  ~isempty(CC) 

frange=linspace(CC(  1),CC(2),  1 0); 
cnt  =  1 ; 

disp(  ’  F requency(Hz)  Pressure(Pa)  ’ ) 


APPENDIX  D.  MATLAB™  FILES 


for  f=frange; 

pp(cnt)=press(p,rtol,NBlock,f,  sigma,  rhop,cs,ca,rhoa,za,zs,z,zO,r,rO,  theta,  thetaO); 

disp([f,abs(pp(cnt))]) 

cnt  =  cnt+1; 

end 

plot(frange,abs(pp),  ’g’) 
hold  on 

xlabel(’Frequency  (Hz)’) 
ylabel(’ Pressure  (Pa)’) 

title(’Soil  Resonance  Predictions  -  Frequency  Variation’) 
end 

D.2.4  In.m 

function  II=In(ka,ks,rhoa,rhos,za,ca,zs,cs, sigma, rhop,f,muv) 
sqeek=csqrt(kaA2-muv); 
queek=csqrt(ksA2-muv); 

II=l/2Hsrhos*(-sin(sqeekHsza).*cos(sqeekH!za)+sqeekH!za)./(sqeek.*cos(sqeek!|sza)  A2)+... 
l/2Hsrhoa*(-queekH!zs+queek*zs.*cos(sqeek*za).A2-cos(queek*zs).*sin(queek*zs)+... 
cos(queek*zs).*sin(queek*zs).*cos(sqeek!|!za).A2 
)./  (cos(queek*zs).A2.*cos(sqeek*za).A2.*queek); 

D.2.5  press.m 

function  pl=press(p,rtol,NBlock,f, sigma, rhop,cs,ca,rhoa,za,zs,z,zO,r,rO, theta, thetaO) 
omega=2*pi*f; 

ks=sqrt((omegaA2+i*omega*sigma/rhop)/csA2); 

ka=omega/ca; 

rhos=  1 1 3  0 ;  %rhop+i  *  sigma/ omega; 

N=0; 

pl=0; 

continN  =  true; 
while  continN 
n=(N+l):(NBlock+N); 

muv  =  mu(n,ka,rhoa,rhos,za,ca,zs,cs, sigma, rhop,f); 
if  muv(l)==kaA2 
muv=muv(2 :  end) ; 
end 

II=In(ka,ks,rhoa,rhos,za,ca,zs,cs,  sigma,  rhop,f,  muv); 

pla=-iV(4Hr).*besselh(0,l,(csqrt(muv)*r0)).*besselj(0,(csqrt(muv).*r)) 

*rho(zO,rhos,rhoa) 

.*Zn(zO,ka,muv,za,ks,zs).*Zn(z,ka,muv,za,ks,zs); 
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pl=pl+sum(pla); 

continM  =  abs(sum(pla)/pl)  >  rtol; 
m=l ; 

while  continM 

p2a=-i./(2*II).*besselh(m,l,(csqrt(muv)*rO)).*besselj(m,(csqrt(muv).5|sr)) 

.  *  cos(m*  (theta-thetaO)) 

*rho(zO,rhos,rhoa).*Zn(zO,ka,muv,za,ks,zs).*Zn(z,ka,muv,za,ks,zs); 

pl=pl+sum(p2a); 

m=m+l; 

continM  =  abs(sum(p2a)/pl)  >  rtol; 
end 

N=N+NBlock; 

continN  =  abs(sum(pla)/pl)  >  rtol; 
end 

D.2.6  propforzl.m 

za=500;  %rough  estimate  of  the  upper  limit  of  the  problem,  units  of  m 
ca=340;  %  speed  of  sound,  m/s 

zs=-0.075;  %rough  estimate  of  lower  limit  of  problem,  units  of  m 
cs=160;  %Attenborough  et  al.  Attached  derivation  from  Buchanan  units  m/s 
cp=270; 
sigma  =  0; 

rhoa=1.2;  %Wikipedia,  density  of  dry  air  at  std  ambient  press  and  temp 

%Density  of  Air  Article,  units  of  kg/mA3 

r0=2; 

z0=2; 

theta0=0; 

NBlock  =  250; 

%rrange=linspaee(0.0 1,1 0,50);  % 
r=l; 

%thetar=linspace(0,90,100); 

theta=0; 

p=[]; 

rtol  =  5e-5; 

rhop=8.1827e7/(cpA2*.27); 

f=319;% 

%frange=linspace(50,300,300); 

z=0; 

cnt  =  1 ; 

zrange=linspace(- 1 ,5, 150); 

%for  r=rrange 
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%for  f=frange; 
for  z=zrange; 

%for  theta=thetar 
omega=2*pi*f; 

ks=sqrt((omega  A2+i*omega*sigma/rhop)/cp  A2); 

ka=omega/ca; 

rhos=1400; 

%pp(cnt)=press(p,rtol,NBlock,f,  sigma,  rhop,cs,ca,rhoa,za,zs,z,zO,r,rO,  theta,  thetaO); 
w(cnt)= 

velpl(p,  rtol,NBlock,rhos,ca,cs,  sigma,  f,  omega,  rhoa,rhop,za,zs,z,zO,r,rO,  theta,  thetaO); 
disp([z,abs(vv(cnt))]);%,abs(vv(cnt))]) 

%disp([z,abs(pp(cnt))]) 

cnt  =  cnt+1; 

end 

plot(zrange,abs(vv),’r’) 

%plot(rrange,abs(vv),  ’g’) 

%Figure 

%plot(frange,abs(pp),’b’) 

%plot(thetar,abs(pp),’b’) 

D.2.7  secant.m 

function  [x,fval]=secant(f,xO,xl  ,tol) 
u=feval(f,xO); 
v=feval(f,xl); 
err=abs(xl-xO); 
while  (err>tol); 
if  (v-u)~=0 

x=x  1  -  v  *  (x  1  -x0)/(  v-u) ; 

xO=xl; 

u=v; 

xl=x; 

v=feval(f,xl); 

err=abs(xl-xO); 

else 

x=xl-(xl-x0)/2; 

xO=xl; 

u=v; 

xl=x; 

v=feval(f,xl); 

err=abs(xl-xO); 

end 
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end 

fval=abs(v); 

D.2.8  velpl.m 

function  vl= 

velpl(p,rtol,NBlock,rhos,ca,cs,  sigma,  f,  omega,  rhoa,rhop,za,zs,z,zO,r,rO,  theta,  thetaO) 
omega=2*pi*f; 

ks=sqrt((omega  A2+i*omega*sigma/rhop)/270.A2); 

ka=omega/ca; 

rhos=rhop+i*sigma/omega; 

N=0; 

pl=0; 

continN  =  true; 
while  continN 
n=(N+l):(NBlock+N); 

muv  =  mu(n,ka,rhoa,rhos,za,ca,zs,cs, sigma, rhop,f); 
if  muv(l)==kaA2 
muv=muv(2 :  end) ; 
end 

II=In(ka,ks,rhoa,rhos,za,ca,zs,cs, sigma, rhop,f, muv); 
if  r<rO 

pla=-i./(4*II).*besselh(0,l,(csqrt(muv)*r0)).*besselj(0,(csqrt(muv)*r)) 

*rho(zO,rhos,rhoa) 

.*Zn(zO,ka,muv,za,ks,zs).*dZn(z,ka,muv,za,ks,zs); 

else 

pla=-iV(4Mr).*besselh(0,l,(csqrt(muv)*r)).*besselj(0,(csqrt(muv)*r0)) 

*rho(zO,rhos,rhoa) 

.*Zn(zO,ka,muv,za,ks,zs).*dZn(z,ka,muv,za,ks,zs); 

end 

pl=pl+sum(pla); 

continM  =  abs(sum(pla)/pl)  >  rtol; 
m=l; 

while  continM 
if  r<rO 

p2a=-i./(2*II).*besselh(m,l,(csqrt(muv)*rO)).*besselj(m,(csqrt(muv)*r)) 

.*cos(mHs(theta-thetaO))*rho(zO,rhos,rhoa).*Zn(zO,ka,muv,za,ks,zs) 

,*dZn(z,ka,muv,za,ks,zs); 

else 

p2a=-i./(2*H).*besselh(m,l,(csqrt(muv)*r)).*besselj(m,(csqrt(muv)*rO)) 

.*cos(m*(theta-thetaO))*rho(zO,rhos,rhoa).*Zn(zO,ka,muv,za,ks,zs) 

,*dZn(z,ka,muv,za,ks,zs); 
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end 

pl=pl+sum(p2a); 
m=m+ 1 ; 

continM  =  abs(sum(p2a)/pl)  >  rtol; 
end 

N=N+NBlock; 

continN  =  abs(sum(pla)/pl)  >  rtol; 
end 

%p=[p,  abs(pl)]; 

vl=-i./(omega.*rho(zO,rhos,rhoa))*pl; 

D.2.9  Zn.m 

function  Z=Zn(z,ka,muv,za,ks,zs) 
kas=csqrt(kaA2-muv); 
kss=csqrt(ksA2-muv); 
if  z>=0 

Z=-(sin(kas*za)./cos(kas*za)).*cos((kas).*z)+sin((kas).*z); 

else 

Z=-sin(kas*za).*sin(kss*zs)/(cos(kas!|!za).*cos(kss*zs)).*(cos(kss*zs)./sin(kss*zs) 

.  *  cos((kss) .  *z)+sin((kss) .  *z)); 

end 


D.3  Eigmovie 

D.3.1  dZn.m 

function  dZ=dZn(z,ka,muv,za,ks,zs) 

sqs  =  csqrt(ksA2-muv);  sqa  =  csqrt(kaA2-muv); 
if  z  >=  0 

%Z  =  cos(sqs*zs).*sin(sqa*(za-z)); 
dZ  =  -sqa.*cos(sqs*zs).*cos(sqa*(za-z)); 
else 

%Z  =  sin(sqa*za).*cos(sqs*(z-zs)); 
dZ  =  -sqs.*sin(sqa*za).*sin(sqs*(z-zs)); 
end 

D.3. 2  eigmovie.m 

clear  all 

zs  = -0.175; 
rhoa  =  1.168; 
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cpa  =  330; 
za  =  500; 
cps  = 160; 
cs  =  75; 

Por  =  0.25; 
rhosn  =  1400; 
nu  =  0.2; 

rhop  =  1.6406e+003; 
sigma  =  lei; 
a  =  0.07; 

T  =  5e2;  %2e6; 
betam  =  T/le3;  %le-2; 
rfreq  =  250; 

rhom  =  T/aA2/(2*pi*rfreq)A2*(2.405)A2; 
rO  =  2;%0.2; 
zO  =  2;  thetaO  =  0; 
z  =  0;  theta  =  0; 

%r  =  linspace(0,l,40); 
r  =1;%  0.0; 

%pts  =  21;  cent  =  248;  freq  =  linspace(cent-5,cent+5,pts); 
pts  =  11;  cent  =  248;  freq  =  linspace(cent-2,cent+2,pts); 
colset  =  [’k’,’b  YcYgYmYr’]; 

%  rO  =  2;  zO  =  2;  thetaO  =  0; 

%  r  =  1 ;  z  =  0;  theta  =  0; 

N  =  5000;  mi  =  0; 

for  n  =  1  :length(freq) 

omega  =  2*pi*freq(n); 

beta  =  sigma/omega/rhop; 

ka  =  omega/cpa; 

ks  =  o  m  ega/cp  s  *  sq  rt(  1  +  i  *  beta ) ; 

rhos  =  rhop*(l+beta*i); 

[mu,kapm]=mumm(N,ka,rhoa,rhos,ks,za,zs); 
plotran  =  300:650;%1 : 1000;%600:800;  % 

%  Figure(2) 

%  subplot(21 1);  plot(real(kapm(plotran))-(2*(plotran)-l)*pi/2,colset(l+mod(n,6))); 
hold  on 

%  subplot(212);  plot(imag(kapm(plotran)),colset(l+mod(n,6)));  hold  on 
%  drawnow 
%  Figure(3) 

%  subplot(21 1); 

%  plot(abs(dZnbnr(z,mu(plotran),za,zs,rhoa,rhos,ka,ks)),colset(l+mod(n,6))); 

%  hold  on 
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%  %plot(l,abs(Znbnr(z,mu(l),za,zs,rhoa,rhos,ka,ks) 
.*Znbnr(zO,mu(l),za,zs,rhoa,rhos,ka,ks)),[colset(l+mod(n, 6)), ’*’]);  hold  on 
%  terms  = 

%  besselh(mi,csqrt(mu(plotran))*rO).*dZnbnr(z,mu(plotran),za,zs,rhoa,rhos,ka,ks) 

.*Znbnr(zO,mu(plotran),za,zs,rhoa,rhos,ka,ks); 

vtemis 

=i/omega/rhoa.*i/2/epsl(mi)./Phib(mu(plotran),za,zs,rhoa,rhos,ka,ks) 

.*besselh(mi,csqrt(mu(plotran))*rO).Hsbesselj(mi,csqrt(mu(plotran))*r)*... 

cos(mi*theta)*rho(z,rhos,rhoa).*Zn(zO,ka,mu(plotran),za,ks,zs) 

,*dZn(z,ka,mu(plotran),za,ks,zs); 

subplot(211);  plot(plotran,real(vterms),colset(l+mod(n,6)));  hold  on 
%subplot(212);  plot(plotran,imag(vtenns),colset(l+mod(n,6)));  hold  on 
drawnow 
%  Figure(4) 

%  plot(abs(besselh(mi,csqrt(mu)*r0).*besselj(mi,csqrt(mu)*r)),colset(l+mod(n,6))); 
hold  on 

fprintf(’f  =  %gHz,  kappa_l  =  %g  +  %gi\n’,freq(n),real(kapm(l)),imag(kapm(l))) 
%p  =  zeros(l,10); 

%  for  m  =  0:9 

%  p(m+l)  =  -sum(i/2/nuO(m)./Phib(mu,za,zs,rhoa,rhos,ka,ks) 
.*besselh(m,csqrt(mu)*rO).*besselj(m,csqrt(mu)*r)*... 

%  sigmaO(zO,rhoa,rhos).*Znb(z,mu,za,zs,ka,ks).*Znb(zO,mu,za,zs,ka,ks) 

*  cos(m*  (theta-thetaO))); 

%  %fprintf(’m  =  %g:  pi  =  %g+  %gi\n’,m,real(p(m+l)),imag(p(m+l))) 

%  end 

%  fprintf(’|p|  =  %g\n’,abs(sum(p))) 

fprintf(’real  sum  =  %g,  imag  sum  =  %g  v  =  %g\n’,sum(real(vterms)), 

sum(imag(vterms)),abs(sum(vterms))) 

pause 

%  Figure(3) 

%  plot(abs(Znbnr(z,mu,za,zs,rhoa,rhos,ka,ks) 
.*Znbnr(zO,mu,za,zs,rhoa,rhos,ka,ks)),colset(n));  hold  on 
end 

D.3.3  epsl.m 

function  Z=epsl(m) 
if  m==0 
Z=2; 
else 
Z=l; 
end 
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D.3.4  Phib.m 

function  Ph  =  Phib(mu,za,zs,rhoa,rhos,ka,ks) 
sqs  =  csqrt(ksA2-mu);  sqa  =  csqrt(kaA2-mu); 

Ph 

=(l/2)*rhos*cos(sqs*zs)  A2.Hs(-cos(sqa*za).*sin(sqa*za)+sqa!|!za)./sqa 
-(l/2)*rhoa*sin(sqa*za).A2.Hs(cos(sqs*zs).*sin(sqs*zs)+sqs*zs)./sqs; 

D.3.5  secant.m 

function  [x,fval]=secant(f,xO,xl  ,tol) 
u=feval(f,xO); 
v=feval(f,xl); 
err=abs(xl-xO); 
while  (err>tol); 
x=x  1  -v.  *  (x  1  -x0)./(v-u); 
x0=xl; 
u=v; 
xl=x; 

v=feval(f,xl); 

err=abs(xl-xO); 

end 

fval=abs(v); 

D.3.6  Zn.m 

function  Z=Zn(z,ka,muv,za,ks,zs) 

sqs  =  csqrt(ksA2-muv);  sqa  =  csqrt(kaA2-muv); 
if  z  >=  0 

Z  =  cos(sqs*zs).*sin(sqa*(za-z)); 
else 

Z  =  sin(sqa*za).*cos(sqs*(z-zs)); 
end 


D.4  Membrane  Problem 

D.4.1  alph.m 

function  TT=alph(zs,ks,muv); 
TT=zs*csqrt(ksA2-muv); 
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D.4.2  Amn.m 

function  y=Arnn(m,ka,ks,muv,zetav,za,zs,rhoa,rhos,a,zO,rO) 
Imnlv=Imnl(ka,ks,muv,zetav,za,zs,rhoa,rhos); 
kk=besselh(m,  1  ,csqrt(muv)*a); 
ll=besselj  (m+ 1  ,csqrt(muv)*a); 
kkl=besselh(m+l,l,csqrt(muv)*a); 

11 1  =besselj  (m,csqrt(muv)*a); 
bessl  =  kk.*ll-kkl.*lll;  %bottom  line 

bess2  =  (-csqrt(muv).*kk  1  ).’*bcssclj(m,csqrt(zctav)*a)+(kk).’*(csqrt(zctav) 
.*besselj(m+l,csqrt(zetav)*a));  %top  line 
%sure  want  to  use  right  divide?  yes 

y=(bess2 .  *Imnlv)\(-i/(2  *  epsl(m)) .  *  csqrt(muv) .  *  (rho(zO,rhos,rhoa) 
*Zn(zO,ka,rnuv,za,ks,zs,rhoa,rhos).*besselh(m,l,csqrt(muv).*rO).*bessl)).’; 

D.4.3  char4z.m 

function  Z=char4z(ka,za,tau,ks,a,k,m,T,rhoa,omega,rhos,zs) 
0/ochar4z(ka,za,kapr,ks,a,k,m,T,rhoa,omega,rhos,zs); 
sqg=kaA2-(tau/za).A2;  %  this  is  zeta 
llamb  =  csqrt(ksA2-sqg); 
sb  1  =besselj  (m,k*  a); 
sb2=besselj  (m+ 1  ,csqrt(sqg)  *a); 
sb3=besselj(m,csqrt(sqg)*a); 

%running  into  problems  here,  too  slow,  line  13,  bessel  function 
IH=a*(sqg4cA2).A(-l).*(csqrt(sqg).*sbl.*sb24c.*besselj(m+l,kH!a).!|!sb3); 

II2=aA2/2  *  (sb3 .  A2+sb2 .  A2)-a./csqrt(sqg) .  *(m.  *  sb2 .  *  sb3) ; 
Z=-rhos*omegaA2./THsl./((kA2-sqg).*sbl).*(cos(tau).*cos(llamb*zs)-rhos.*tau 
./(za*rhoa*llarnb).*sin(tau).*sin(llamb*zs)).*(besselj(m,a*csqrt(sqg)).H!IIl-sbl.*II2)+... 
II2.*(-llamb.H!cos(tau).*sin(llambHszs)-rhos.*tau./(rhoa*za).*sin(tau).!|!cos(zs*llamb)); 

D.4.4  dZn.m 

function  dZ=dZn(z,ka,muv,za,ks,zs) 

sqs  =  csqrt(ksA2-muv);  sqa  =  csqrt(kaA2-muv); 
if  z  >=  0 

%Z  =  cos(sqs*zs).*sin(sqa*(za-z)); 
dZ  =  -sqa.*cos(sqs*zs).*cos(sqa*(za-z)); 
else 

%Z  =  sin(sqa*za).*cos(sqs*(z-zs)); 
dZ  =  -sqs.*sin(sqa*za).*sin(sqs*(z-zs)); 
end 
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D.4.5  epsl.m 

function  Z=epsl(m) 
if  m==0 
Z=2; 
else 
Z=l; 
end 

D.4.6  Imnl.m 

function  Y=  Imnl(ka,ks,muv,zetav,za,zs,rhoa,rhos) 

M  =  length(zetav);  N  =  length(muv); 
xiv=za*csqrt(kaA2-muv).’*ones(l,M);  %1 
iotav=zs*ones(N,l)*csqrt(ksA2-zetav);  %n 
alphv  =zs*csqrt(ksA2-muv).’  *ones(  1  ,M); 
ttau=za*ones(N,l)*csqrt(kaA2-zetav); 

%  Yl=  -za*rhos*(ttau.*cos(ttau).*sin(xiv)-xiv.*cos(xiv).*sin(ttau)) 
./(cos(ttau).*cos(xiv).*(ttau.A2-xiv.A2)); 

% 

%  Y112  =  (zs*sin(xiv)./(za*cos(ttau).*cos(xiv).*cos(alphv) 

.*(alphv.A2-iotav.A2))); 

% 

%  Y121  =  zs*rhos*ttau.*(cos(alphv).*cos(ttau)-cos(ttau).*cos(iotav)); 

%  Y122  =  rhoa*za*(alphv.*sin(ttau).*sin(alphv)-iotav.*sin(iotav).*sin(ttau)); 

% 

% 

zsan=l/2*(-2*sin(xiv).A2*rhoa*zs.*(alphv+cos(alphv).*sin(alphv))./(cos(xiv).A2 

.*cos(alphv).A2.*alphv)).A(l/2)+l/2*2A(l/2)*(rhos*za*(-sin(xiv).*cos(xiv)+xiv) 

./ (xiv.  *  cos(xiv) .  A2)) .  A(  1/2); 

% 

% 

phin=l/2*2A(l/2)*(rhos*za*(-cos(ttau).*sin(ttau)+ttau)./(ttau.*cos(ttau).A2)).A(l/2) 
+  l/2*(-(-2*rhosA2*ttau.A2.*cos(iotav*zs).*sin(iotav*zs) 
+2*rhosA2*ttau.A2.*iotav*zs 

+4*rhos*ttau.*tan(ttau).*cos(iotav*zs).A2*rhoa*za.*iotav 

+2Hstan(ttau).A2HsrhoaA2*zaA2.*iotav.A2.*cos(iotav*zs).*sin(iotav*zs) 

+2*tan(ttau).A2*rhoaA2*zaA2.*iotav.A3*zs-4*rhos*ttau.*tan(ttau)*rhoa*za.*iotav) 

,/(rhoa*iotav.A3)).A(l/2)/za; 

% 

%  Y12  =  Y121  +Y122; 

% 

%  Y2=Y12.*Y112; 
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% 

%  Y=(Yl+Y2)./(zsan.*phin); 

%new  IMNL 

Y  =  zs*sin(xiv).*(-zs*rhos.*ttau.*cos(iotav).*cos(ttau)-iotav.*rhoa.*sin(iotav) 
.*sin(ttau)*za+zs.*rhos.Hsttau.*cos(alphv).*cos(ttau)+rhoa*alphv.*sin(alphv) 
.*sin(ttau)*za)./(za*(-iotav.A2+alphv.A2))+zaHscos(alphv).Hsrhos.5|s(-ttau.*cos(ttau) 
.*sin(xiv)+xiv.*sin(ttau).*cos(xiv))./(-ttau.A2+xiv.A2); 

D.4.7  iota.m 

function  00=iota(zs,ks,zeta); 

00=zs*csqrt(ksA2-zeta); 

D.4.8  numb.m 

function  nO=numb(za,ka,ks,zs) 

n0=round(za.*sqrt(ka.A2-ks.A2+(pi/(2*zs)).A2)/pi)+1000; 

D.4.9  phip.m 

function  y=phip(ka,zetav,z,za,rhoa,rhos,zs,ks) 
ttaua=csqrt(kaA2-zetav); 

%phin=l/2*2A(l/2)*(rhos*za*(-cos(ttau).*sin(ttau)+ttau)./(ttau.*cos(ttau).A2)).A(l/2) 
+  l/2*(-(-2*rhosA2*ttau.A2.Hscos(iotav*zs).*sin(iotav*zs)+2*rhosA2*ttau.A2.*iotav*zs 
+4*rhos*ttau.*tan(ttau).*cos(iotav*zs).A2*rhoa*za.*iotav 
+2  *  tan(ttau) .  A2  *  rhoaA2  *  zaA2 .  *  iotav  A2 .  *  cos(iotav  *  zs) .  *  sin(  iotav  *  zs) 
+2Hstan(ttau).A2*rhoaA2*zaA2.*iotav  A3*zs-4*rhos*ttau.H!tan(ttau)*rhoa*za.5|siotav) 
,/(rhoa*iotav.A3)).A(  l/2)/za; 

%y=ttau/za.*(cos(ttau).H!cos(ttau/za*z)+sin(ttau).*sin(ttau/za*z))/phin; 
y=-c  sqrt(kaA2  -zetav) .  *  sin(ttaua  *  (z+za)); 

D.4.10  propforzl.m 

za=500;  %rough  estimate  of  the  upper  limit  of  the  problem,  units  of  m 
ca=330;  %  speed  of  sound,  m/s 

zs=-.075;  %rough  estimate  of  lower  limit  of  problem,  units  of  m 
cs=170;  %Compressional  Wave  Speed 
cp=  75;  %Shear  Wave  Speed 
sigma  =  0; 

rhoa=  1.168;  %Wikipedia,  density  of  dry  air  at  std  ambient  press  and  temp 

%Density  of  Air  Article,  units  of  kg/mA3 

r0=2; 
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z0=2; 

theta0=0; 

%rrange=linspac  e(0 . 0 1 , 1 0 , 1 00) ;  % 

r=0.0; 

theta=0; 

p=[]; 

rtol  =  5e-5; 

rhop=  1 453  ;%8 . 1 827e7/(cpA2  *  .27); 
frange=linspace(80,300, 150); 
rhom=6 1 84;%5980.04; 
a=.07; 

T=2e6; 

betam=1000; 

%zrange=linspace(  -  .075, .075, 100);% 
z=0; 

NBlock=250; 

N  =  numb(za,ka,ks,zs); 

NN=1:N; 

MM=0; 
cnt  =  1 ; 

%  omega=2*pi*f; 

%  ks=sqrt((omega.A2+i*omega*sigma/rhop)./csA2); 

%  ka=omega/ca; 

%  k=csqrt((rhom  *  omega .  A2+i  *betam  *  omega)/T) ; 

%  rhos=rhop+i*sigma./omega; 

%  muv  =  mu(N,rhos,rhoa,ks,ka,za,zs); 
%II=Imnl(ka,ks,muv,zetav,za,zs,rhoa,rhos); 

YY1=[]; 

%for  z=zrange 
%for  r=rrange 
for  f=frange 
omega=2*pi*f; 

ks=sqrt((omega  A2+i* omega* sigma/rhop) ./csA2); 
ka=omega/ca; 

k=csqrt((rhom*omega  .^-K^betam^omegayT); 
rhos=1400; 

muv  =  mu(N,rhos,rhoa,ks,ka,za,zs); 

YY=velol(omega,rhoa,muv,r,MM,ka,ks,N,k,T,rhom,za,zs,rhos,a,zO,rO,theta,z); 

disp([f,abs(YY)]); 

YY 1 = [YY 1  ,(YY)] ; 
end 

plot(frange,abs(YY  l),’r’) 
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%plot(rrange,abs(YYl),’b’) 

%plot(zrange,abs(YYl),’g’) 

%y=Amn(M,ka,ks,muv,zetav,za,zs,rhoa,rhos,a,zO,rO); 

%zrange=linspace(0,5,200); 

%  for  f=frange; 

%  %for  z=zrange; 

%  pp(cnt)=velo  l(p,rtol,N,f, sigma, rhop,cs,ca, cm, rhoa,za,zs,z,zO,r,rO, theta, thetaO, a, T); 
%  disp([f,abs(pp(cnt))]) 

%  %disp([z,abs(pp(cnt))]) 

%  cnt  =  cnt+1; 

%  end 

%  %plot(zrange,abs(pp),’-m’) 

%  plot(frange,abs(pp),  ’r’) 

D.4.11  secant.m 

function  [x,fval]=secant(f,xO,x  1  ,tol) 
u=feval(f,xO); 
v=feval(f,xl); 
err=abs(xl-xO); 
while  (err>tol); 
x=x  1  -v.  *  (x  1  -x0)./(v-u); 
x0=xl; 
u=v; 
xl=x; 

v=feval(f,xl); 

err=abs(xl-xO); 

end 

fval=abs(v); 

D.4.12  velol.m 

function  v  1  =velo  1  (omega, rhoa,muv,r,m,ka,ks,N,k,T,rhom,za,zs,rhos, a, z0,r0, theta, z) 
vl=0; 
for  ii=m 

zetav=zeta(N,ka,za,ks,a,k,m,T,rhoa,omega,rhos,zs,rhom); 

vla=l./(i*omegaH!rhoa)*besselj(ii,(csqrt(zetav)5|!r)) 

*Amn(ii,ka,ks,muv,zetav,za,zs,rhoa,rhos,a,zO,rO)H!cos(ii*theta) 

*phip(ka,zetav,z,za,rhoa,rhos,zs,ks); 

vl=vl+sum(vla); 

end 
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D.4.13  xi.m 

function  UU=xi(za,ka,muv); 

UU=za*csqrt(kaA2-muv); 

D.4.14  zeta.m 

function  [zeta,kap]=zeta(N,ka,za,ks,a,k,m,T,rhoa,omega,rhos,zs,rhom) 
ftntol=le0; 
sectol=5e-4; 

%  width=pi/8; 

%  Nlook=min(N,floor(  1 12  *  ( 1 +2  *ka*za/pi))+ 100); 
lasteig=0 ;  count=0 ; 
k  1 =zeros(  1  ,N)  ;k2=zeros(  1  ,N) ; 
kapr=linspace(i*le-3,500*i,2000); 

y=imag(char4z(ka,za,kapr,ks,a,k,m,T,rhoa,omega,rhos,zs)); 

if  y(l)>0 

l=find(y<0); 

else 

I=find(y>0); 

end 

if  -iscmptyO) 
count=count+l; 
k  1  (count)=kapr(I(  1 )- 1 ); 
k2(count)=kapr(I(  1)); 
end 

while  count  <=  N 

kapr=linspace(lasteig+ 1  e-3  ,lasteig+2  *pi,200) ; 

y=real(char4z(ka,za,kapr,ks,a,k,m,T,rhoa,omega,rhos,zs)); 

if  y(l)>0 

I=find(y<0); 

else 

I=find(y>0); 

end 

if  (m==0)|(((kapr(I(l)-l)  <  za*ka)  &  (kapr(I(l))  >  ka*za))) 

count=count+l; 

k  1  (count)=kapr(I(  1 )- 1 ); 

k2(count)=kapr(I(  1)); 

lasteig=k2(count); 

else 

lasteig=za*ka+ 1  e-3 ; 

end 

end 
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%  if  count  <  N 

%  kl((count+l):N)  =  kappa((count+l):N)-width; 

%  k2((count+l):N)  =  kappa((count+l):N)+width; 

%  end 

[kap,fv]=secant(@(kap  1)  char4z(ka,za,kap  1  ,ks, a,k,m,T,rhoa, omega, rhos, zs),k  1  ,k2,sectol); 
%  if  N  >  Nlook 
%  NN  =  (Nlook+l):N; 

%  II  =  find(abs(kap-kappa(NN))  >  ftntol); 

%  for  nn  =  II 
%  ii  =  nn+Nlook; 

%  [kapc,xfv]=fmins  1 04( ’prob  1  obj f ’ ,[real(kappa(ii)),imag(kappa(ii))], [] ,pi/4 
/real(kappa(ii)),ka,ks,k,m,za,zs,rhoa,rhop,sigma,a,T,rhom); 

%  kap(ii)  =  kapc(l)+i*kapc(2); 

%  fv(ii)  =  xfv; 

%  fprintf(’%g  ’,ii) 

%  if  mod(ii,10)  ==  0 
%  fprintf(’\n’) 

%  end 
%  end 

%  fprintf(’\n’) 

%  end 

zeta=kaA2  -(kap/ za) .  A2 ; 

D.4.15  Zn.m 

function  Z=Zn(z,ka,muv,za,ks,zs,rhoa,rhos) 

sqs  =  csqrt(ksA2-muv);  sqa  =  csqrt(kaA2-muv); 

normz=l;%-rhosH!cos(sqs.*zs).!|!(cos(sqa.*za)-l)./sqa-rhoa*sin(sqa.H!za).*sin(sqs.*zs)./sqs; 
if  z  >=  0 

Z  =  cos(sqs*zs).*sin(sqa*(za-z))./normz; 
else 

Z  =  sin(sqa*za).*cos(sqs*(z-zs))./normz; 
end 
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This  section  includes  only  variables  which  are  used  in  more  than  one  section. 

A,  B,C, ...  -  arbitrary  constants  (may  or  may  not  have  subscripts). 

a  -  diameter  of  the  membrane. 

c  -  speed  of  sound  through  a  medium. 

ca  -  propagation  speed  through  air. 

cs  -  propagation  speed  through  porous  soil. 

D  -  region  occupied  by  the  membrane  (Ch.  4  and  5). 

Dm  -  weight  function  for  the  non-homogeneous  membrane  problem  (Ch.  4). 
F  -  radial  function  in  cylindrical  coordinates  for  separation  of  variables. 

/  -  frequency. 

G  -  angular  function  in  cylindrical  coordinates  for  separation  of  variables. 

Q  -  Green’s  function.  (Ch.  4) 

g  -  gravitational  constant,  g  =  9.8  m/  s2  =  32.2  ft  /  s2. 

H  -  humidity.  (Ch.  2) 

Hm  -  Hankel  Function. 

7 H  -  porosity.  (Ch.l) 
i  -  imaginary  number,  i  =  \/—l. 

J  -  flux.  (Ch.  1) 

Jrn  -  Bessel  function  of  the  first  kind  of  order  m. 
k  -  wave  number,  k  —  - 

C 

/C  -  structure  constant.  (Ch.  1) 
ka  -  wave  number  in  air.  ka  =  — 

Ca 

ks  -  wave  number  in  soil.  ks  = 

m  -  index  value,  m  —  0, 1,  2, ...  (unless  otherwise  indexed) 
n  -  index  value,  n  —  0, 1, 2, ...  (unless  otherwise  indexed) 

P  -  pressure.  (Ch.  1) 

P  -  total  pressure,  P  =  P  +  p.  (Ch.  1) 

p  -  small  perturbation  in  the  pressure,  acoustic  pressure.  (Ch.  1). 
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p  -  pressure  function.  (Ch.  2). 
pa  -  pressure  function  in  the  air. 
pe  -  external  pressure.  (Ch.  5) 

Pinc  -  incident  pressure  field.  (Ch.  5) 
ps  -  pressure  function  in  the  soil. 

Pscat  -  scattered  pressure  field.  (Ch.  5) 

Q  -  rate  at  which  density  is  created  per  unit  volume.  (Ch.l) 

R  -  small  element  of  volume.  (Ch.  1) 
r  -  radial  component  in  cylindrical  coordinates 
r0  -  radial  component  of  the  source. 

S  -  entropy.  (Ch.  1) 

S  -  total  entropy,  S  =  S  +  s.  (Ch.  1) 
s  -  small  perturbation  in  entropy.  (Ch.  1) 

T  -  temperature.  (Ch.  1). 

T  -  total  temperature,  T  =  T  +  r.  (Ch.  1) 

T  -  time  function.  (Ch.  2). 

T  -  membrane  tension  (Ch.  4  and  5). 
t  -  time. 

u  -  fluid  flow  velocity.  (Ch.  1) 
u  -  vertical  displacement  of  the  membrane.  (Ch.  5) 

V  -  a  small  volume.  (Ch.  1) 
w  -  velocity  of  flow.  (Ch.  1) 

x  -  the  horizontal  coordinate  in  the  standard  xyz  axis  configuration. 

Ym  -  Bessel  function  of  the  second  kind  of  order  m. 
y  -  the  vertical  coordinate  in  the  standard  xyz  axis  configuration. 

Z  -  height  function  in  cylindrical  coordinates  for  separation  of  variables. 

Za  -  height  function  in  the  atmosphere  (0  <  z  <  za). 

Zs  -  height  function  in  the  soil  (zs  <  z  <  0). 

z  -  the  height  coordinate  in  the  standard  xyz  axis  configuration. 

z  -  the  height  component  in  cylindrical  coordinates. 

z0  -  the  height  component  in  cylindrical  coordinates  of  the  source 

za  -  the  height  limit  for  the  model  of  the  atmosphere  (maximum  atmospheric  altitude). 

z,s  -  the  depth  (negative  height)  for  the  rigid  substrate. 

X  -  the  weight  function  for  the  depth  equation,  equation  (4.19). 

5  -  small  perturbation  in  the  density.  (Ch.  1). 

8  -  Dirac  delta  "function". 

£  -  small  value  added  to  the  radial  component  in  the  Green’s  Function. 
em  -  values  in  Mathews  and  Walker’s  Green’s  Function.  (Sec.  4.3,  Ch.  5) 

7  -  ratio  of  specific  heat  at  constant  pressure  to  specific  heat  at  constant  volume,  T 
ks  -  compressibility  of  the  elastic  membrane  (Ch.l). 
k  -  substitution  variable,  n  =  za\Jks  —  //. 

A  -  arbitrary  constant  for  separation  of  variables,  found  to  be  A  =  n2,  where  n  =  0, 1,  2, ... 
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fi  -  arbitrary  constant  for  separation  of  variables,  found  to  be  /i  =  k2 

Qa  -  substitution  variable  for  \Jk2  —  (  .  (Ch.  5) 

0.5  -  substitution  variable  for  ^Jk2  —  (.  (Ch.  5) 

(Ir„  -  substitution  variable  for  f  x  (z)  Z2  (z)  dz. 


2 


2s 

(j)  -  depth  condition  above  the  membrane. 

d)a  -  depth  condition  above  the  membrane  in  the  atmosphere. 

(j)s  -  depth  condition  above  the  membrane  in  the  soil. 
p  -  density. 

p  -  total  density,  p  =  p  +  5.  (Ch.  1) 
pa  -  air  density. 

pp  -  effective  density  of  the  porous  soil. 
ps  -  soil  density. 

7 r  -  pi,  3.14159...,  the  ratio  of  the  circumference  to  the  diameter  of  the  circle. 
C  -  substitution  variable.  (  =  fir2. 

(  -  eigenvalues.  (Ch.  5) 
a  -  flow  resistivity. 

9  -  angular  component  in  cylindrical  coordinates. 

60  -  angular  component  in  cylindrical  coordinates  of  the  source, 
r  -  small  perturbation  in  temperature.  (Ch.l) 
r  -  substitution  variable,  r  =  za  \/k2  —  Qnn.  (Ch.  5) 
uj  -  angular  velocity,  u  =  2nf 


